Reports of the Department of Geodetic Science 

Re port No. 245 

JNASA-CE- 156^970) DESIGN OF EXPESIHENT FOE N78-22440 

EASTS ROTATION AND BASELINE PARA2ETER 
DETERMINATION FROM VERY LONG BASELINE 

INTEREER OHETRY {Ohio State Uni-v. Research Dnclas 

.Foundation) 123 p BC A06/HF API • CSCL 08E G3/43 18679 J 

DESIGN OF EXPERIMENT FOR EARTH ROTATION 
AND BASELINE PARAMETER DETERMINATION 
FROM VERY LONG BASELINE INTERFEROMETRY 



National 


by 

Athanasios Dermanis 


Prepared for the 

Aeronautics and Space Administration 
Washington, D.C. 


Grant No. NGR 36-008-204 
OSURF Project No. 783820 



The Ohio State University 
Research Foundation 
Columbus, Ohio 43212 


February, 1977 



Reports of the Department of Geodetic Science 
Report No. 245 


DESIGN OF EXPERIMENT FOR EARTH ROTATION AND 
BASELINE PARAMETER DETERMINATION FROM 
VERY LONG BASELINE INTERFEROMETRY 


by 

Athanasios Dermanis 


Prepared for the 

National Aeronautics and Space Administration 
Washington, D.C. 

Grant No. NGR 36-008-204 
OSURF Project No 3820-A1 


The Ohio State University 
Research Foundation 
Columbus, Ohio 43212 


February, 1977 



PREFACE 


OSURF Project 3820-A1 is under the supervision of Professor Ivan I. 
Mueller, Department of Geodetic Science, The Ohio State University, and under 
the technical direction of Mr. James P. Murphy, Special Programs, Office of 
Applications, Code ES, NASA Headquarters, Washington, D.C. The contract is 
administered by the Office of University Affairs, NASA, Washington, D.C. 20546. 

The VLBI data and other pertinent information received from Mr. Peter F. 
Mac Doran of the Jet Propulsion Laboratory in connection with this investigation is 
gratefully acknowledged. 


n 



TABLE OF CONTENTS 

Page 

Preface ... . ii 

1. Introduction 1 

2. VLBI— Basic Technique and Mathematical Model 3 

3. Earth Potation Parametrization 6 

4. Linearization of Observation Equations 11 

5. Coordinate System Definition and Inner Constraints 16 

6. Remarks Related to Design of Experiment Concepts 19 

7. Simulation and Adjustment Philosophy 25 

8 . Arrangement of Experiments 28 

9. Results and Conclusions 35 

References 67 

Appendix A: Inner Constraints and Their Relation to the "Geometry” 

of the Pseudoinverse of an Operator 68 

Appendix B: Algorithm for First-Order Partitioned Linear Regression 

Including Inner Constraints 72 

Appendix C: Computer Programs 76 



LIST OF FIGURES „ 

Pag 

2.1 The geometry of VLBI observations 4 

6.1 Geometric loci of parameter sensitivity vectors 23 

9.1 Geographic locations of stations and baselines in 

Experiments 1 through 10 34 

9.2 Recovery of earth rotation parameters, Experiments 1-10 39 

9.3 Recovery of baseline lengths and angles. Experiments 1-10 50 

9.4 Variation of earth rotation parameter recovery with steps ize and 

the total time interval of observations (Experiment 2) 61 

9.5 Variation of baseline length and angle recovery with steps ize and 

the total time interval of observations (Experiment 2) 62 

9.6 Variation of earth rotation parameter recovery with earth step 63 

9.7 Variation of baseline length and angle recovery with earth step 65 

A.l The geometry of the pseudoinverse operator 70 

LIST OF TABLES 

6.1 Sensitivity of Baseline - Radio Source Extreme Configurations 

with Respect to Various Parameters 21 

8.1 Present and Prospective VLBI Station Locations and Their 

Approximate Coordinates 29 

8.2 Fictitious Stations Used in Experiments 6 and 9 30 

8.3 Experiments and Participating Baselines 31 


IV 



1 . Introduction 


Very Long Baseline Interferometry (VLBI) is one of the new techniques 
which will probably dominate geodesy and geophysics in the near future . Its main 
advantage lies in the fact that it brings the accuracy of direction measurements to 
a level previously possible only for range measurements . This closes the gap 
between powerful range determination techniques such as laser ranging and the 
much less accurate determination of directions through photographic tracking of 
artificial earth satellites. 

The technique is geometric in the sense that the observations ar e independent 
of the gravity field of the earth. However, the ’’orbits" of the observed extragalactic 
radio sources with respect to an earth-fixed system are dominated and perturbed by 
the rotation of the earth with respect to inertial frame . This allows the determination 
of polar motion, precession-nutation and length-of-th e-day variations, and the 
technique becomes also "dynamic” in this respect. 

The capability of determining the geometry of a network of stations within 
a short time interval and with a centimeter level accuracy also allows the, study 
of the variation of network geometry with time caused by earth tides and other 
periodic or secular station drifts . 

The obvious importance of VLBI for both geodetic and geophysical applications 
is only limited by the effect of atmospheric refraction and-unwanted noise associated 
with instrumentation [Shapiro and Knight, 1970, Section 3]. At present the develop- 
ment of the technique is at the stage of establishing its capabilities by direct com- 
parison with classical surveys in relatively short baselines where geophysical 
effects are negligible, and of making the system more transportable. 

In the present work the effect of nonwhite instrumentation noise, atmospheric 
refraction and earth tides are ignored, and the emphasis is on the possibility of 
recovering earth rotation and network geometry (baseline) parameters. It is 
assumed that nonwhite instrumentation noise is either eliminated or separately 
determined through calibration except for a linear effect which cannot be 
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distinguished from clock errors included in the adjustment parameters. Errors 
in the determination of atmospheric refraction are assumed to be either absent 
or uncorrelated in which case their effect is incorporated in the variance of the 
considered instrumentation -white noise. The effect of' earth tides on the variation 
of station coordinates is assumed to be calculated from separate information and 
subtracted from the actual observations. 

The numerical simulated experiments performed here are setup in an 
environment where station coordinates vary with respect to inertial space according 
to a simulated earth rotation model similar to the actual but unknown rotation of 
the earth. 

In Chapter 2 the basic technique of VLBI and its mathematical model are 
presented. The parametnzation of earth rotation chosen is described in Chapter 3, 
and the resulting model is linearized in Chapter 4. Cartesian station coordinates 
are not estimable quantities in the analysis of VLBI observations. It is possible to 
use a model where only estimable quantities are considered (see, for example, 
[Arnold, 1974]) . However, the choice of Cartesian coordinates leads to a simpler 
model, and the problem of coordinate system definition is resolved in Chapter 5 
with the use of inner constraints. A simple analysis of the geometry of the observa- 
tions in Chapter 6 leads to some useful hints on achieving maximum sensitivity of 
the observations with respect to the parameters considered. The basic philosophy 
for the simulation of data and their analysis through standard least squares 
adjustment techniques is presented in Chapter 7. 

The main objective of the present work is the exploration of the capabilities 
of VLBI for the recovery of earth rotation and baseline parameters . For this 
purpose, a number of characteristic network designs based on present and candidate 
station locations is chosen in Chapter 8. The results of the simulations for each 
design are presented in Chapter 9 together with a summary of the conclusions . 
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2. VLBI— Basic Technique and Mathematical Model 


Very Long Baseline Interferometry (VLBI) is a technique for observing 
time intervals and/or their derivatives at two widely separated antennas between 
the two instants of reception of the same wavefront emitted from extragalactic 
radio sources. The fundamental difference with respect to classical radio (short 
baseline) interferometry is in the fact that the ante nna s are not directly connected 
by cables for the comparison of the received signals . Sufficiently stable atomic 
frequency standards at each site make it possible to record the signals as functions 
of time. The recorded signals are then crosscorrelated afterwards m a computer. 
The fact that maximum correlation occurs when the two recordings are shifted in 
time by an amount equal to the wave travel time between the two instances of 
reception, makes it possible in principle to determine this time interval as well 
as its time derivative. Extragalactic radio sources are sufficiently far away for the 
wavefront to be considered as a plane and thus from the known velocity of the wave 
(velocity of light) the geometric distance corresponding to the observed time delay 
can be calculated. 

The crosscorrelation procedure is described in detail in a series of 
articles by Thomas (1972a, 1972b and 1973), and we shall be concerned here only 
with the geometry of the observations. In Fig. 1, XYZ is a geocentric Cartesian 
reference frame fixed with respect to the radio sources (assumed to be an inertial 
frame). A certain wavefront arrives at station i at epoch ti, when the position 
vector of the station is Xi(ti), and at the station j at a slightly different epoch t 2 
when the position vector of this station is X j(t^ . The true time delay t 2 - ti is the 
travel time corresponding to the projection Dij of the retarded baseline X j(ts) - Xi(ti) 
on the direction of the radio source. If e is the unit vector in the radio source 
direction and c the velocity of light, 

Dij = c (t a - fe) = [Xj(to) - Xi(fa)l * e (2-1) 
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Since (tg - t x ) is a very short time interval, the time variation of the position vector 
Xj(t) may be assumed to be solely due to the earth's rotation. Therefore 


X,(fe) = Xj(ta) + 



ti 


t 2 

O(t) x Xj(t) dt 


( 2 - 2 ) 


where fi(t) is the rotational velocity vector of the earth. £2(t) is practically constant 
over the short time interval (ti - tg), and thus, by setting Q(t) = fi(ti) 

Xj(ts) = Xj(ti) + (t 2 -ti) [0(60 X Xj(t!)] (2-3) 

Setting Bij(t) = Xj(t) - Xi(t) and Vj(t) = 0(t) X Xj(t), one obtains 

D U = Bij(ti) * e - - Djj V-j(ti) ■ e (2-4) 


If approximate values of 0(tx), Xj(ti) and e are known with sufficient accuracy such 
that the corresponding uncertainty in the retardation effect 

6D tJ = - Du Vj * e (2-5) 

Q 

is negligible, one obtains, by setting dij = Dij + 6 and switching to matrix 
(column vector) notation, for an observation from the ij baseline to the radio source 
p at epoch t*, the much simpler model 

dijpk = Bij (1*) e p (2-6) 

Assume that d iJlJt has been corrected for effects giving rise to discrepancies between 
modeled and true distance, e.g., atmospheric errors , relativistic effects," aberration, 
except for clock errors. Assuming that the two local oscillators at stations i and j 
differ only by a constant offset C 0 ij (at some epoch to) and a linear drift Cij, the 
model becomes 

dijpk = Blj(t k ) e p + e[Coij + Cu(t k - t 0 )] (2-7) 

Since the inner product of two vectors is invariant with respect to the choice of 
reference frame, the model may be rewritten 

dijpk = b T u e p (tk) + c [C 0 ij + Cij(tk-t 0 )] (2-8) 
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with 

blj = [Xj-xj, y s -y lf Zj-zj] (2-9) 

where bjj and e p are the baseline and the radio source direction unit vector, 
respectively, with respect to an earth-fixed reference frame x,y, z. One also has 

e* (%) = M(i*) e p (2-10) 

where M(t) is the earth rotation matrix of transformation from the inertial to the 
earth- fixed frame. The model with respect to the earth-fixed reference frame 
becomes 

= blj M(t k ) e p + c[Cotj + C u (tk - to)] (2-11) 

Further development of the model depends on the particular parametrization chosen 
for the rotation of the earth. 


3. Earth Rotation Parametrization 

The most straightforward way of parametrizing earth rotation is that of 
expressing the inertial to earth-fixed system transformation matrix M(t) in terms 
of three rotation angles at every epoch t. The most logical choice seems to be that 
of the Euler ian angles <p, 0, jii, defined as 

M(t) = Ba [<p(t)] Ba. [0(t)J Bs [0( t)] (3-1) 

where Ri, R s , Rs denote conventional rotation matrices about the x, y,z axes, 
respectively. The rotating earth can be viewed as a dynamical system whose 
state evolves in time according to a second-order differential equation of the 
form 

= f(E, t) (3-2) 

where E = [<p, 0, 0] T . If such equations of motion could be exactly known and solved, 
their solution 

E(t) = F(E (to) , E(t 0 ), t) ( E = dE/dt) (3-3) 


6 



would provide a parametrization of the transformation matrix 


M(t) - Mftfb, 6o, 0o, <Po, do, 0O to, t) 
in terms of the six initial values (integration constants): 


<Po — ^5(to) 0o — ®(to) 0O ” 0(to) 
• = * 2 . fl. = d§ = M 


= 


dt j to 


dt | to 


dt to 


However, such an approach is not practical because of uncertainties surrounding 
the current knowledge of the earth's rotation. The alternative is the representation 
of the functions <p (t), 0(t), 0(t) in terms of a finite number of parameters to be 
estimated from the observations. Such representations as polynomials, trigonometric 
series, etc. can sufficiently approximate a wide class of functions, provided that 
appropriate number of terms is included The efficiency of the approximation can 
be judged only when the function to be approximated is known over some time interval, 
while in our case an unknown function is to be approximated. The unknown approxima- 
tion error (representation error, modeling error) may well be significant compared 
to the observational errors, thus affecting the validity of the final parameter error 
estimates based upon real data analysis. 

Among various finite parameter representations, the simplest perhaps is 
that of representing a continuous bounded function over an interval by a step function, 
which has a constant value over each subinterval, mto which the original time 
interval has been partitioned. 

To formally represent a step function, we introduce the characteristic or 
indicator function Xa of a set A, defined as follows: 

/ 1 if xEA 

Xa (x) - < (3-5) 

( 0 if x<£A 

Let [a,b] be the original interval and 

a = to < t x < t a < ... < Vx < t* “ b 
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If x i denotes the characteristic function of the set tj] (i = 1, 2, . . . , n-1) and 

Xn that of [tn_i, tn] , a step function with respect to the above partition of [a, b] can 
be written in the form 

n 

y(t) “i Xi (t) (3-6) 

A continuous bounded function f(t) can be approximated arbitrarily well by a step 
function y(t) with an appropriate choice and increase in the number of subintervals. 
In the case of this investigation, the functions to be approximated are the Eulerian 
angles cp(t), 9(t), 0(t), and the effectiveness of step function approximation maybe 
increased if the known parts of <p, 0, 0 are filtered out. Let <p°(t) , 0°(t), 0°(t) be the 
"reference" Eulerian angles as computable from classical astronomy. Then the 
earth rotation transformation matrix may be represented as 

M(t) - Ik[<p 0 + 6<p] + 60] R 3 [0°+60] (3-7) 

where 6co (t), 6 0(t) and 50 (t) are step functions . An alternative parametrization 
may be of the form 

M(t) = M°(t) 5M(t) = R 3 (p°) Ri(0 p ) R 3 (0°) R 3 (6<P) Ei(60) R 3 (50) (3-8) 

Traditionally, earth rotation representations are separated into three parts 
with the help of the instantaneous rotation vector co(t) related to the Eulerian angles 
through Euler's geometric equations 


00 x 

= sin 0 sintp 

0 

+ 

cos <p 

0 


COy 

= sin 0 cosip 

0 

+ 

sin <p 

0 

(3-9) 


= cos 0 

0 

+ 


4- < p 



where a) x , co y , are the rotation vector components with respect to the earth- 
fixed system. The variation in the direction of the vector co with respect to the 
inertial system constitutes precession and nutation, while polar motion is the 
corresponding variation with respect to the earth-fixed system. The third part of 
earth rotation is the variation in the length of the vector co, i.e., the variation of 
the rotational velocity (length of the day) £2 = || do|[ . 
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To obtain a parametrization of earth, rotation where precession-nutation, 
polar motion and rotational velocity are separated, we introduce an intermediate 
moving geocentric Cartesian reference frame x, y', z . The z-axis direction 
coincides at any epoch with the direction of the instantaneous rotation axis, while 
the x-axis direction is arbitrary. The transformation from the X, Y, Z (inertial) 
to the x, y', z! frame may be expressed in terms of three rotations: 

Rs (© 2 ) Ra(S) Rj (H) 

and the transformation from x, y', z' to the x,y, z (earth -fixed) frame can be 
similarly expressed as 

Ri(-T?) R 2 (- t) Rs(©i) 

Setting © = 0i + © 2 , the total transformation from the inertial to the earth-fixed 
frame becomes 

M = Ri(-7?) Re(-£) Rs(©) Re(S) Ri(H) ( 3-10 ) 

where 17 , £, ©, H, H are all functions of time, (Note that the ©1 and -© 2 components 
of © cannot be separated, and this justifies the arbitrariness of the x-axis direction. ) 
In order to introduce the instantaneous rotational velocity £2 in the model, noting that 
the R 3 (©) rotation is about the instantaneous rotation axis (diurnal rotation), one 
may set 

t 

0 (t) = ©o +J £ 2 (T) dT ©o = ©(to) (3-11) 

to 

The three angles H, H, 0 O represent the traditional precession-nutation, and the 
angles 7) and £ the polar motion. The angle (© - © 0 ) is similar to GAST (Greenwich 
Apparent Sidereal Time) . The functions 7 ), £, £2, H, H can now be approximated by- - 
step functions. To improve the approximation, we shall include in the transformation 
M the traditional precession-nutation based on current knowledge, as follows: 

M - Ri(-7?) Rj>(~£) Rs(©) Rs(H) Ri(ll) R° (3-12) 

where the known part of precession and nutation may be computed from [Mueller, 

1969, Ch. 41. 
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R° = Ri(-e-Ae) R 3 (-A0) Ri(e) R 3 (-z) Re(0) Rs (-Co) 


(3 13) 


In the above expression, z, 0, Co are the Newconib components of precession, and 
A ij) and Ae the nutations in longitude and obliquity by-Woolard. The transformation 
R° defines the traditional "true” frame of reference, an intermediate, noninertial, 
moving reference frame X', y', Z f of known orientation, such that the Z -axis is 
at any epoch near the instantaneous rotation axis. The angles E, H, ©o are thus 
small corrections to the reference precession-nutation parameters used in (3-13), 
and as such they may be considered as the transformation parameters between the' 
"traditional (Newcomb/Woolard) true" frame and the "actual true" frame. 

There are two alternatives in the step function representation of the angle © . 
Gne may represent O as a step function with constant value Stover the i^ 1 time 
interval (step), in which case . ' . 

„ .i-i - 

©(t) = ©o+X) (t* - tk-o + a t (t - t t -a) (3 _ 14) 

k 1 ' ti_i < t £ tj 

A second choice is to approximate,© over each step by a function of the form 

©i(t) = ©oi + £2 t (t-ti_ x ) t,_i < t< t t (3-15) 

This is a weaker representation from the adjustment point of view because of the 
larger number of parameters (© 0 i, i = 1, 2, . . . , N vs. 0 O ), but it leads to estimates 
of variances of the estimates of the step function constants 7h, 0 o t, £2 t , E t , Ht, 
which are uniform (equal) over all steps (tj-i, ti). On the contrary, the first 
representation leads to variances which are smaller for steps .in the- middle and 
larger for steps at the ends of the original time span of the available observations. 
Over each step the model for the VLBI observations becomes 

dij* = blj Ri(-t?) Bb K) Rs [©o + a (t* - t,)] Rs (E) Ri (H) R° e P + 

+ c [C otJ + Cu(t k - tb)] (3-16) 

where 7), 0 O , £2, E, H are constants over the considered step (the subscripts 

identifying the steps have been dropped for simplicity), t k is the epoch of observation, 
and to is the beginning of the step containing t k . 
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4. Linearization of Observation Equations 


For the time span of available observations considered in this study 
(maximum of one month), the effect of precession-nutation on the experimental 
design (relative radio sources - stations configuration) is negligible. One can 
therefore drop the known transformation matrix R° from the model without any 
significant change on the simulation results. This simplifies somewhat the equations 
following, but it must be understood that appropriate modifications must be made 
when analyzing real data by including the effect of R°. 

The angles £, 17 , E, H are very small since both the z and Z axes (actually 
the noninertial Z -axis on account of ignoring R°) are or can be chosen to be near 
the rotation axis. The usual approximations, cos a = 1, sin a = a, a 2 = 0, for a 
small angle a may be used, leading to 

/ 1 '° 

Ri (~V) RaK) = 10 1 

\~t V 

and 

l 1 0 

R 2 (E)R x (H) -(01 
\ H -H 

The model now becomes 

dij* = ( r J “ r i ) T S T <£,??) Ra (©o + ft Tfc) S(E, H) e P + 

+ c [Coij + C tJ Tfc] (4-1) 

where T k = t* - tp; e p - [cos 6 P cos 0! p , cos 6 p sinarp, sin 6 P ] T , 0! p and6pare 
the right-ascension and declination of the p^ 1 radio source with respect to the 
X, Y, Z inertial frame; r t = [x l5 yi, zi] ' , xi, yi , zi are the coordinates of the 
i^ 1 station with respect to the x,y,z earth-fixed frame; and 

c 

S(a, b) •= I 0 
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The available observations dt Jjk form a vector of observations Lb, corresponding to 
a vector of observables L a = F(X a ), where X a is a vector of unknown parameters . The 
unknowns are the station coordinates x t , y t , z 1} i = 1, 2, . . ., N s , the earth rotation 
parameters £, 77 , 0 O , Q, E, H (one set of six per each step), radio source coordinates 
a*, 6 P , p = 1, 2, . . . , Nrs, and clock synchronization parameters C 0 ij, Cjj, i, j = 1, 2, 

. . , N s . 

The dfjpfc components of the vector La are linear with respect to most of the 
parameters with the exception of © 0 , O, Op, and 5 P . A completely linearized model 
can be obtained by means of approximate values X° of the parameters X a . Expansion 
in Taylor series and neglect of second- and higher -order terms leads to 


L a 


F(X°) + 


d E(X a ) 
BXa 



(X a - X°> 


(4-3) 


The actual observations L b differ from the observables L a because of the presence of 
unknown observations errors V, according to L a - L b + V. Setting A = (SF/oX a ) | x o, 

L° - F(X°), X =X a - X°, and L = L° - L b , the linearized observation equations are 
obtained: 

V = A X + L (4-4) 

A and L are known, while X and V are unknown. The vector V is assumed to be an 
outcome of a random vector with E[V] = 0 and known covariance matrix E[V V T ] = S. 

Application of standard least squares techniques leads to a solution X, 
minimizing the quadratic form V T PV where P = S 1 , which satisfies the "normal 
equations, " 

(A t PA)X = -A t P L or NX = U - (4-5) 

To compute the elements of the desing matrix A, the partial derivatives of dj jpt with 
respect to the parameters are required. Introducing the notation, - © Q + Cl T k , 

Xkp =#k - ttp, =Xj - x t , yjj = yj - y u z i} == Zj - zj, and dropping subscripts where 
confusion is not likely to arise, the partial derivatives are as follows: 
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3d_ 

3xj 

3d 

Syj 

3d 

3 Zj 


3d 

3xi 

3d 

3yi 

3d 

5zj 


cos 6 cos X + sin 6 (- E cos 0 + H sin 0 + 4) 

- cos 5 sin X + sm 6 (E sin 0 + H cos ip - rj) 

sin 6 - cos 6 (4 cos X + V sin X ~ S cos O' + H sin a) 


= x tJ [sin 6 + cos 6 (E cos or - H sin a - 4 cos X)] + 

+ Yij [17 cos 6 cos X] + 

+ Zij [- cos 6 cos X - sin 6 (- E cos 0 + H sin ip + 4)1 

— = y 13 [- sin 6 + cos 6 (4 cos x + sin X _ cos a + H sin a)] + 

+ Zij [- cos 6 sin x + sin 6 (E sin0 + H cos ip - 77)] 


= xjj [- cos 6 sin X + sin 6 (E sin!/) + H cos ip)] + 

+ yij [~ cos 6 cos X + sin 5 (E cos ip - H sin !/>)] + 
+ Zij [cos 6 (4 sin X ~ V cos X)] 

3d 3d 3d 3d 

30 “ Tk 3©o * 3C 0 ij " C ’ 3Cij ° Tk 


= xjj j-cos ip cos 5 + cos 6 [4 cos 0! + cos ip (- E cos a + H sin a)] j + 

+ yij | sin!/) sin 6 + cos 6 [- 7] cos a + sin 0 (E cos a - H sin O')] j + 
+ Zi) [cos 6 cos a + sin 6 (4 cos 0 + T) sin 0 - E)] 


= xij [sin 6 sin 0 + cos 6 sin Of' (E cos 0 - H sm 0 - 4) ] + 

+ yij [sin 6 cos 0 + cos 6 sin a (- E sin 0 - H cos 0 + r?) ] + 
+ Zij [- cos 6 sin 0! + sin 6 (- 4 sin 0 + 1? cos 0 - H) ] 

dd 

— = Xu [ cos 6 sin X] + yu [cos 6 cos x] + 

+ Zij cos 6 (- 4 sm X + T) cos X - E sin oc - H cos a) 
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OQ * 

gg = xij [~ sin 6 cos X + c °s 6 (- H cos ip + H sin$ + £ ) ] + 

+ yij [sin 6 sin x + cos 6 (H sin $ + H cos ij) - rj ) ] + 

+ Zij [cos 6 + sin 6 (£ cos X + V sin X - E cos or + H sin a) ] 

(4-6) 


A second type of possible observations is the time derivative of d tjl)c 

fijkp = = SI (r i- rO T S t (|,77) P 3 B3{0 o + SI r k ) S(H,H) e p + c C„ 

(4-7) 

where 


P 3 



1 

0 

0 



(4-8) 


The related partial derivatives of fi j!cp with respect to the parameters are as follows: 


J_ Sf 
SI d£ 


_1 3f 
Sr] 


1_ df 

o as 


[£ cos 6 sin x] - yij [V cos 5 sinX] 

+ zij [cos 6 sin x - sin 6 (E sin 0 + H cos $) ] 

yij [cos 6 (- £ sin x + V cos X) ] 

+ z 1} [- cos 6 cos X + sin 6 (H cos ip - E sin ip) ] 


xij [sin 6 sin ip - cos 6 sin ip (- S cos or + H sin or) ] 

+ yij [sin 6 cos ip - cos 6 cos ip (- H cos a + H sin or) ] 
+ zij [sin 6 (- £ sin0 + r\ cos ip) ] 


jY = x tJ [cos ip sin 6 - cos 6 sin or (S sin ip + H cos ip) ] 

+ yij [- sin ip sin 5 + cos 6 sin a (- E cos ip + H sm ip) ] 
+ /-ij [sm 6 (- £ cos ip - 7) sin ip )] 
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i a f 


Cl d© 0 


Xij [- cos 6 cos X + sin 6 (E cos 0 - H sin ip) ] 

+ Yu [ eos 5 sin X - sm 6 (- H sin if) - H cos if)) ] 
+ zij [cos 6 (£ cos X + V sin X) ] 


1 

a f 

1 

a f 

n 

a Xj 

n 

dxi 

i 

a f 

i 

a f 

n 

3 yj 

" Cl 

dyi 

i 

aj_ 

i 

a f 

n 

azj 

n 

azi 

i 

Cl 

aj 

aa 

= xij [cos 6 i 


■pr -r — = - cos 6 sin x + sin 6 (H sin if) + H cos 0) 


= - cos 6 cos X + sin 6 (S cos if) - H sin 0) 


= cos 6 (£ sin X - V cos X) 


+ zij [cos 6 (- £ cos x - V sin x) ] 


_1 cif 

a as 


Xij [sin 6 sin X + cos 6 (E sin 0 + II cos 0) ] 

+ yij [sm 6 cos X + cos 6 (E cos 0 ~ H sin 0) ] 
+ zij [sin 6 (- £ smx +7) cos X) ] 


M - r il 
an " k a© 0 ’ 


ii- = o — 

ac o£J ’ acij 


(4-9) 
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5. Coordinate System Definition and Inner Constraints 


The combination of all observations leads to the linearized observation 
equations 

*Vi = n Au + n Li (5-1) 

Minimization of V T PV leads to a solution satisfying the normal equations [Uotila, 1967] 
NX = U _ (5-2) 

where 

N = A t PA and U = -A T PL. 

Because of the lack of coordinate system definition, N is singular with rank 

(N) = r < u, where u is the number of parameters and s = u - r is the rank deficiency 

of N. 

Originally s = 9, corresponding to the six degrees of freedom in the definition 
of the earth-fixed coordinate system (three for origin position and three for orientation) 
plus the three degrees of freedom in the inertial system definition (orientation only) 

Since a radio source catalogue of a certain accuracy is assumed to be available, 
radio source coordinates are treated as observations (weighted unknowns) rather than 
as unknown parameters. The inertial frame is therefore defined through the catalogued 
radio source coordinates and the remaining rank deficiency of N becomes s = 6. A 
unique solution to the normal equations may be obtained if in addition a set of minimal 
linear constraints are imposed on the parameters, of the form 

C T X = 0 (5-3) 

where C is a u X s matrix and rank (C) = s. 

Among the various possible solutions to the normal equations, the unique one 
given by X - N* - U, where N + is the pseudoinverse of N, has the following properties 
[Blaha, 1971]: X T X = min. and trace N + = min. In view of the interpretation of N + 
as the variance -covariance matrix of the parameters, the second property makes 
the solution optimal in the sense that smaller variances of the parameters provide a 
better representation of parameter related estimable quantities in terms of the 
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non estimable parameters. To avoid the use of pseudo inverse computation algorithm, 
one may resort to a particular set of minimal constraints called inner constraints 
[Blaha, 1971], 

E t X = 0 ' (5-4) 

which leads to the same solution for X as the one obtained with the use of the 
pseudoinverse N + . It can be shown that such a set of inner constraints can be 
obtained with the help of a u X s matrix E (rank (E) = s) satisfying 

N E - 0 (5-5) 

An algebraic type of proof can be found in [Pope, 1971] which settles the truth of the 
matter but throws little light on the interrelation between inner constraints and 
pseudo inverse. A different type of proof is given in Appendix A of this work, based 
on the geometry of the operator represented by the matrix N. 

In view of the fact that AE = 0 => NE = A T PAE = 0, a set of inner constraints 
has been analytically constructed using six independent solutions to the set of equations 

A J y = 0 j = 1, 2, . . u (5-6) 

where A j denotes the j^ 1 row of A. We assume that the approximate values of 
7), £, H, H are zero in the computation of partials of observations with respect to 
the unknowns, thus obtaining simpler analytical expressions for the elements of A. 
Letting the order of the unknowns in X to be the following 

X T = [dr?, d£, d© 0 , dtt, dxi, dyi, dzi, ..., dx N , dy N , dz N ] 


one iLicty 


' y ~ [wi, w 3 , ws, w 4 , a a , &, y-i, y u .... 

corresponding to a di jkp observation, one has 




For a row of A 


9d 

M 


9d 


dd 

ad 

a ^d 

ad 

wi y— + 


+ w 3 

+ 

o 

© 

<V 

w 4 


+ - + 
dxi 

A ^ + 



Sd 

+ ft 

Sd 


ad 




+ 

a, — 

3xj 

3yj + 

y j 

azj 

= 0 




(5-7) 
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With the analytical expressions of the partials and after a considerable algebraic 
effort, one arrives at six independent solutions yi, 1 = 1, 2, . . . , 6, and the inner 
constraint matrix E is formulated by setting 

E = [yi y 2 ... y 6 ] 

The final result is 


d 7] dr} d©o dfl dx x dyi dzi 

0 0 0 0 1 0 0 

0 0 0 0 0 1 0 

0 0 0 0 0 0 1 

1 0 0 0 0 -z° y° 

0 1 0 0 z? 0 -x? 

0 0-10 -y? x° 0 


dxu dy N dz N 

10 0 
0 10 

0 0 1 

A 0 0 

o -Z N y N 

o _ o 

z N 0 -Xfi 

-yN xn o 

(5-8) 


where x°, y°, z° are the approximate values of station coordinates used in the 
conputation of the design matrix A More explicitly there are two sets of constraints. 
The first one 



(5-9) 


"defines" the origin of the earth- fixed system, while the second one 
/-dr? \ / 0 -z? y? \ /dxi \ 

( "d£ j = S j z° 0 -xi j I dyi J (5-10) 

\d© 0 / \-y? x? 0/ \dz t / 


"defines" its orientation. 

i 

In the discussion above, the problem of epoch and time scale definition for the 
set of station clocks has been ignored The parameters C 0 ij and Ctj refer only to 
relative offsets and drifts between the clocks at stations i and j. 
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Usually epoch and time scale are provided by the readings of one of the 
clocks (master clock), and they are transformed to other station clocks through the 
parameters C QD j, C n j where the master clock is at the station. 

Another way to provide epoch and timescale is the use of inner constraints 
defining a fictitious master clock. We may shift to a new set of parameters by 
setting 

Cotj = C 0R j - C 0 pi Cij = Cpj - Cpi (5-11) 

The new parameters are offsets C 0 pj and drifts Cpiwith respect to the fictitious 
master clock, and the additional inner constraints are 

dC oP i = 0 and ^ dC p i = 0 (5-12) 

i i 



6. Remarks Related to Design of Experiment Concepts 

Since the usual objective of a geodetic experiment is the estimation of 
certain parameters, an experimental design may be considered optimal when the 
errors in parameter estimates are minimized in some certain sense . Within the 
linear least squares model, the a posteriori variance -covariance matrix of parameter 
estimates is given in general by 

S* - Ob 8 (A t PA) + = ct 8 Q + (6-1) 

+ 

where Q denotes the pseudoinverse of a square matrix Q, P is the weight matrix 
and Oo the a posteriori estimate of the variance of unit weight. The ’’errors” in 
parameter estimates therefore depend on the observational accuracy and on the 
design matrix A. The design matrix A = A(X°) itself depends on the approximate 
values of parameters X°. The determination of optimal X° values minimizing some 
risk function (e g. , trace Sx) is usually referred to as the "configuration problem, " 
or the "first-order optimal design" [Grafarend, 1974] On the assumption that the 
risk function chosen is insensitive to small changes of X°, the approximate values 
of parameters may be identified with the time ones. In a similar way the small 
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angles 77 , £, 5, H in our model may be considered to be zero. Of the remaining 
parameters 0 O depends on the choice of reference frames, Q is approximated by 
the fixed average rotational velocity of the earth, while in the case of transcontinental 
baselines the choice of radio source coordinates a-p, 6 p .is strongly restricted by 
observability conditions. Therefore, the first-order optimal design is practically 
limited to the determination of the optimal configuration of the network of observing 
stations . 

The order of magnitude of elements in a row of the design matrix A determines 
the sensitivity of the corresponding observation with respect to each of the parameters . 


Next a short geometric interpretation of these partials is given to provide possible 
"hints" towards the determination of the optimal design. 

The partial derivatives with respect to 77 , £, © 0 , E, H (earth rotation 
parameters) and ap, 6 P are proportional to the baseline length. Indeed, they all 
contain linear terms in xjj = rij x?j , y^ = rij y°j and Zjj = rj 3 zij, where rij is the 
length of the ij baseline and x°j, y?j, zij are the components of the unit vector r°j 


in the direction of the baseline. It is therefore possible to define "sensitivity per 

1 Sd 

unit of baseline length" by means of — t-x where |3 stands for any of the'earth 

T’lj °p 

rotation or radio source coordinates parameters. It can be shown that 


rij 5/3 


(rV Sp 


( 6 - 2 ) 


i.e. , that the sensitivity per unit of baseline length is the projection of some 
vector Sg (called the parameter ]3 sensitivity vector) on the direction of the baseline. 
It can be easily verified that the parameter sensitivity vectors are given by- 

= - Ri (ir/2) P yz Rs (X) e& 

S £ = - R 2 (77-/2) P xz Rs (X) e& 

S 0 ^ — Tjf - P xy Ra (X +lt/2) e& 

s v = Ra #>) R 2 ( 77 -/ 2 ) P xz R 3 (- a) e 6 ' 

S H = Rs(0) Rl (tt/ 2) P yz Ra(-a) e 6 

S a = ^ Ra (X + 77/2) P xy e 6 

S 6 = Rs(X) R 2 ( 77 / 2 ) e& (6-3) 
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where 

ee = [cos 6, 0, sin 6] T , e p = R 3 (-a) e& p and 

/I 0 0\ /0 0 0\ /I 0 0\ 

Pxy = ( 0 1 0 J , P yz - o 1 0 ) , Pxz = ( 0 0 0 

\0 0 0 / \0 0 1 / \0 0 1 / 

are projection matrices on the xy, yz and xz planes, respectively. The.position 
of the sensitivity vectors varies with time, and their loci, all contained within the 
unit sphere, are depicted in Fig. 6. 1. Maximum sensitivity occurs when the baseline 
is parallel to the corresponding parameter sensitivity vector, while an observation 
is insensitive to a parameter when the baseline and parameter sensitivity vector 
are perpendicular. With this in mind and with the help of Fig. 6. 1, Table 6. 1 is 
constructed, summarizing the sensitivity to various parameters for extreme cases 
of baseline and radio source directions . 


Table 6.1 

Sensitivity of Baseline - Radio Source Extreme Configurations 
with Respect to Various Parameters 


Parameters 

Baseline Parallel to 
Equator (z u = 0) 

Baseline Parallel to Rota- 
tion Axis (xjj = y u = 0) 

Equatorial 
Radio Sources 
(6 = 0°) 

Polar 

Radio Sources 
(6 = 90°) 

Equatorial 
Radio Sources 
(6 = 0°) 

Polar 

Radio Sources 
(6 = 90°) 

V 

no 

yes 

yes 

no 

E, H 

no 

yes 

yes 

no 

©o, a 

yes 

no 

no 

no 

6 

no 

yes 

yes 

no 
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Note that observations in the fourth colu m n are impossible due to observability- 
conditions (a baseline parallel to the rotation axis must be located near the 
equator and observe equatorial radio sources only) . 

The sensitivity with respect to Cl as well as to Cij is proportional to the time 
span r k of the available observations while that to Cij and C 0 ij is independent of the 
baseline - radio source configuration. 

Sensitivity with respect to station coordinates is independent of baseline 
length. If e x , e y , e z denote unit vectors in the directions of the x, y, z axes, 
respectively, the relevant partial derivatives under the same approximations are 


3d 

3xj 


S c , 


3d 

Syj 


By Sc, 


3d 

3zj 


el Sc 


where S c - R 3 (~X) eg. The sensitivity of observations with respect to station coordinates 
is determined by the projection of the coordinate sensitivity vector S c on the coordinate 
axes. The locus of S c is depicted in Fig. 6.1. It can readily be seen that observations 
to equatorial radio sources are insensitive to the z station coordinates, while observa- 
tions to polar radio sources are insensitive tox andy station coordinates. 

All the above remarks are of a rather general nature and simply provide 
,T hints" towards the design of an optimal station configuration for the VLBI experiment, 
which also depends on which of the parameters in particular is to be determined. 
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Fig .6.1 Geometric loci of parameter sensitivity vectors , 
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7. Simulation and Adjustment Philosophy 


The objective of this simulation is the determination of standard deviations of 
parameter estimates relevant to VLBI observations for a number of characteristic 
network designs. In addition to this design (relative configuration of station - radio 
sources and arrangement of observations in time), the a posteriori standard devia- 
tions of parameters also depend on a number of other factors. These include the 
accuracy of observations, their total time span, and the modeling of earth rotation 
and relative station motions in the adjustment. Since step function models are used 
here, the time interval over which earth rotation parameters are considered to be 
constant ("earth step") and the time interval over which station coordinates are con- 
sidered constants ("station step") become also variables. It has also been assumed 
that a radio source coordinate catalogue is available, and its accuracy is one addi- 
tional variable affecting final results. Summarizing, the problem is to examine the 
variation of standard deviations of the following parameters: 

i. station coordinates 

ii. earth rotation parameters 
with respect to the following independent variables: 

a. observational noise 

b. time span of observations 

c. earth step 

d. station step 

e. accuracy of radio source catalogue 

f. station - radio source configuration, arrangement of 
observations in time. 

With respect to the last item above, three different variables appear in the same 
group. This is justified in view of the fact that once a transcontinental network 
design has been decided upon, observability conditions impose strong limitations 
on both the choice of radio sources to be observed and die arrangement of observations 
in time. 
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The fact that observations are not instantaneous, but averages over some 
time interval (integration time of 7-10 minutes), imposes also limitations on the 
maximum density of observations in time. 

The choice of earth and station step in real data analysis depends on the 
frequency content of the corresponding time functions . In view of the aliasing 
effect [Koopmans, 1974], the smaller the step size, the higher the spectral resolution 
of the corresponding functions . On the other hand, larger step size leads to a 
smaller number of parameters to be estimated and a "stronger" adjustment with 
more degrees of freedom . The choice of optimal step size in the analysis of real 
data constitutes a very interesting problem worthy of a separate study. 

In the present simulation study, step size has been considered as a variable, 
while the corresponding functions generated in the simulation of observations are 
continuous rather than step functions . However, precautions have been taken so 
that the total variation T of the function f(t) in question over each step time 
interval I is insignificant compared to observational noise, hi fact, 

V f Ti - max [f(t)] - min [f(t)J ^ a b /N i - 1, 2, ...,n 

’ tei t ten 

where a b is the standard deviation of the observations, and N is a large -integer 
(specifically here N = 10). 

In view of a 7-10-minute integration time, and allowing for antenna motion, 
observations (time delays and their time derivatives) are uniformly spaced at 
15-minute intervals over the total time span considered. Their corresponding 
standard deviations CT b and o/ are taken to be of a constant ratio oyW = 1/2, where 
Ob is in nsec and o/ in nsec/hr. This value is typical, but the ratio may be somewhat 
different depending on antenna and observed signal characteristics. -The choice of a 
fixed ratio reduces the number of considered parameters by one with a corresponding 
significant reduction in computational effort. 

Once a certain observational pattern has been decided upon, a simulation 
program (VLBI SIMULATOR) produces perfect observations using the simplified 
model (Eq. 4-1). This program requires as input station and radio source true 
coordinates and the pattern of the observations (baseline - radio source - epoch of 
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observation) . The station coordinates are taken to be constants, while the values 
of earth rotation parameters at the epoch of observation are provided by a sub- 
routine (SUBROUTINE EARTH). 6b is fixed to a realistic value for the beginning 
of the observations to. The effect of traditional precession-nutation is absent from 
the simulation as well as from the adjustment model for computational simplicity. 
The inertial frame thus is taken to be the true equatorial system of the epoch to so 
that its Z-axis is near the rotational axis. For the same reason, the earth-fixed 
frame is taken to be the CIO/BIH zero-meridian terrestrial system. 

The perfect observations are processed through a program which adds 
"Observational errors" drawn from a Gaussian population with zero mean and 
desired variance. 

For the adjustment of the observations, two separate programs are used. 

The first one (PROGOl) processes observations of total time span equal to both 
earth and station step (one day). An option is available for processing any desired 
number of different radio source catalogue accuracies in the same run. 

The second program (PROG02) processes observations of time span equal 
to the station step only, which in turn is taken to be an integer multiple of the 
earth step. It is, however, possible to process observations using unequal earth 
steps within the total span of available observations. Within each earth step the 
linearized observations equations are of the form 

Vi = AjX + AjXj + Li 

where X refers to parameters common to all earth steps (station and radio source 
coordinates) andXi refers to parameters particular to each earth step (earth 
rotation parameters) . 

The total set of such individual observation equations is combined in a unified 
solution through a "second-order partitioned linear regression" scheme, extended 
to incorporate inner constraints. The detailed algorithm is presented in Appendix 
B. The related programs and subroutines are in Appendix C. 
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8. Arrangement of Experiments 


To obtain more realistic results, station locations considered here are 
not completely arbitrary, but limited to present and prospective locations of 
VLBI stations. Such candidate station locations are presented in Table 8. 1. In 
the same table the stations actually used in this study are indicated by asterisks . 

Two criteria have been considered in the selection. First, stations close to plate 
boundaries, i.e., in geophysically active areas, have been excluded to avoid possible 
comparatively large station drifts deviating from those considered in our simula- 
tions. Secondly, from groups of stations close to each other, only one station has 
been included on the assumption that final results would be similar for a nearby 
station. It is recognized that these criteria probably exclude the most active stations 
currently in operation, but this should not influence the general conclusions sought. 

An experiment consists of a group of stations simultaneously participating 
in the observations. In Table 8. 3 the station groups for the experiments considered 
in this study are presented. After a certain station group has been selected, the 
arrangement of observations is one observations every 15 minutes, which interval 
and choice of the radio sources are limited by observability conditions. The 
observability condition naturally imposed is that the observed radio source must be 
above the horizon at both ends of the observing baseline. To avoid large atmospheric 
refraction errors, the zenith distance of the radio source is limited, in general, to 
be less than z ^ = 45°, with some exceptional cases where z^ = 60°. 

Such observability conditions limit the observable part of the sky to within the 
intersection of two cones with the local station verticals as axes and z teX as the 
vertex angle. For transcontinental baselines, the region of observability is small 
and the choice is between observing a certain quasar or a neighboring one without 
significant effect on the experimental design. 

In setting out a certain sequence of observations, the choice between 
possible alternatives has been guided by the hints obtained in Chapter 6, with 
emphasis on sensitivity with respect to earth motion and station coordinate parameters. 
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Table 8.1 


Present and Prospective VLBI Station Locations and 
Their Approximate Coordinates 


No. 

Location 

Longitude 

Latitude 


Kashima, Japan 

137° 

35° 

1 

Canberra, Australia * 

149 

-35 

2 

Kauai, Hawaii * 

200 

22 

3 

Fairbanks, Alaska * 

212 

65 


Goldstone, California 

243 

35 


Algonquin, Canada 

277 

51 


Greenbarik, Virginia 

282 

38 

4 

Haystack, Massachusetts * 

288 

41 


Santiago, Chile 

289 

-34 


Arecibo, Puerto Rico 

294 

18 

5 

Sao Paulo, Brazil * 

313 

-24 

6 

Madrid, Spain * 

356 

40 


Bonn, Germany 

7 

51 

B 

Onsala, Sweden * 

17 

59 

8 

Johannesburg, So. Africa * 

28 

-26 


Crimea, USSR 

34 

45 


* indicates stations considered in this study 
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Table 8.2 


Fictitious Stations Used.in.Exper-iments 6-and 9' 


No. 

Location 

Longitude 

Latitude 

9 

equatorial 

0° 

0° 

10 

ft 

60 

0 

11 

II 

120 

0 

12 

tt 

180 

0 

13 

It 

240 

0 

14 

IT 

300 

0 

15 

meridian 

0 

0 

16 

n 

0 

60 

17 

tt 

180 

60 

18 

n 

180 

0 

19 

it 

180 

-60 

20 

tt 

0 

-60 
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Table 8.3 


Experiments and Participating Baselines 
(Station Nos . Refer to Numbering in Table 8.1 or 8.2) 


Experiment 

No. 

Participating Baselines 

1 

2-3, 2-4, 3-4 

2 

3-4, 3-7, 4-7 

3 

4-5, 4-6, 5-6 

4 

5-6, 5-8, 6-8 

5 

1-2, 2-3, 2-4, 3-4, 3-7, 4-5, 4-6, 4-7, 5-6, 5-8, 
6-7, 6-8 

6 

9-10, 10-11, 11-12, 13-14, 14-9 

7 

2-4, 4-5 

8 

2-4, 4-6 

9 

15-16, 16-17, 17-18, 18-19, 19-20, 20-15 

10 

3-4, 4-7 
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An effort lias also been made to observe radio sources of a variety of 
declinations to avoid the critical configuration appearing when, for example, all 
radio sources are of the same declination. From a larger set of available radio 
sources, only a rather uniformly distributed set has been considered. The 
criterion for inclusion has been the use of the radio source in previous experi- 
ments so that its appropriateness for VLBI observations is assured. 

Since a radio source catalogue of a certain accuracy is assumed to be avail- 
able, no effort has been made to optimize radio source coordinate recovery. As a 
result, this study is of no "astrometric" value, and no significant improvement of 
radio source coordinate accuracy can be found in the results of the adjustment 
performed. 

One sidereal day can be taken as the "unit" for an observational pattern since 
a sequence of observations performed over one day can be repeated over the next 
days. For observations of a total time span of a number of days, the same daily 
pattern is repeated in the simulation. Various experimental designs are compared 
to each other on the basis of one day of observations where earth rotation and baseline 
parameters are considered as constants over the whole day. 

The choice of experiments performed has been primarily directed towards 
the possibility of recovering earth rotation parameters from a minimum number of 
observing stations . The minimum number is the three stations necessary for the 
definition of an earth -fixed system through minimal or inner constraints. Experi- 
ments 1, 2, 3 and 4 are eases of such minimal "triangle networks" where all the sides 
of the triangle are observing baselines . 

A weaker desing (Experiments 7, 8, 10) is a three station configuration where 
only two of the triangle sides are observing baselines . In particular. Experiment 10 
is a counterpart of Experiment 2, such that the effect of removing one baseline (other 
design parameters remaining the same) can be studied. The opposite to such minimal 
three-station designs is provided in Experiment 5 where all eight stations participate 
in the observations. 
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In addition to these "realistic" station locations, two "fictitious" designs 
are also considered, their common characteristic being the closure of a network 
of stations on a great circle around the earth. Experiment 6 is such a network 
of uniformly spaced stations along the equator, while Experiment 9 is of a 
similar design along a meridian. The coordinates of stations involved in these 
experiments are listed in Table 8.2. 

Fig. 9. 1 depicts the geographic locations of die networks in the above- 
mentioned experiments. 
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Fig. 9.1 Geographic locations of stations and baselines in 
Experiments 1 through 10. 
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9 . Results and Conclusions 


The results of the simulations are presented in three parts. The first part 
(Figs. 9.2- 9. 3) deals with the comparison of the ten designs described in the 
previous chapter, on the basis of one day of observations. Fig. 9.2 depicts the 
recovery (standard deviations) of earth rotation parameters (4, rj, E, © 0 , fl), 
while Fig. 9-3 shows the recovery of network geometry in terms of baseline lengths 
and angles (estimable quantities) . 

Standard deviations a for each parameter are given as a function of the 
standard deviation Oq of radio source coordinates and for various values of obser- 
vational precision CT b . Values used in the simulation are cr b , a q = 1, 5, 10, 20, 50 cm. 

To achieve uniformity, ct q is expressed as arc length on a reference sphere with a 
radius of 6371 km (!"« 30 m). 

The second part (Figs. 9. 4-9 5) is a study of the improvement of the standard 
deviations of the parameters with the total time interval T of available observations. 
Station coordinates are assumed to be constant over T while earth rotation parameters 
are constant only over the sub intervals of duration At (earth step). The ratio p = o7cr b 
is given for each parameter as a function of T for various values of At. In the simu- 
lations, the following values have been used: T = 1, 2, 3, 4, 5, 10, 30 days; At = 6, 

12, 24 hours. Experiment 2 has been chosen for this study. Radio source coordinate 
standard deviations are assumed to be ct.q = o b for all cases. However, the variation 
of the results with < 7 Q is insignificant for T ^ 5 days. 

The third part (Figs .9.6 and 9.7) is a study of the variation of p = cr/CT b 
with At (earth step) based on one day of observations (T = 1 day) for some of the other 
designs (Experiments 1, 3, 4, 6, 9, 10) is shown. 

Since the results of the simulations arc fully presented here, one may draw 
his own conclusions. However, the following general remarks are mild e in an attempt 
to summarize the results: 

(1) Earth Rotation Parameters from One-Day Observations 

For one day of observations and the recovery of earth rotation parameters, 
the following are noted: 

(a) Three-baseline "triangle" designs (Experiments 1, 2, 3, 4 in Fig. 9.2) 
are capable of recovering polar motion (£, rp with a precision of the same order of 

magnitude as that of the obsei vations, even for low radio source accuracies (Oq = 50 cm). 
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For two-baseline designs (Experiments 7, 8, 10 in Fig. 9.2), polar motion recovery 
is of the same quality as for the "triangle" station configurations. However, the 
removal of one baseline in Experiment 10 results in an increase of polar motion 
standard deviations by a factor of 1. 5 compared to Experiment 2. The more complex 
design of Experiment 5 offers no dramatic improvement in polar motion recovery. 

An essential improvement can be seen in the equatorial design of Experiment 6 and 
even more in Experiment 9 which is the "best, " but admittedly not very practical 
design for polar motion recovery. 

(b) The standard deviations of the (precession -nutation) parameters E and 

H are practically the same for all two- or three -baseline designs (Experiments 1, 2, 
3, 4, 7, 8, 10). In particular, the removal of one baseline in Experiment 10 with 
respect to Experiment 2 has no effect on the recovery of E and H. The best 
recovery of E, H can be seen in Experiments 5 and 9. Although the station con- 
figuration in Experiment 6 is of the same shape as in Experiment 9, the change in 
their position (equatorial vs. meridian) has a great effect on E, H standard deviations 
Experiment 6 has the weakest recovery, while Experiment 9 provides the strongest one. 

(c) The standard deviation of © 0 varies with the station configurations . It is 
highest for Experiments 4, 8, 7; medium for Experiments 10, 3, 5; and smallest 
for Experiments 1, 2, 9, and especially 6 (best). 

(d) Recovery of 0 varies strongly even for designs of the same type. It is 
generally weak for two-baseline designs (Experiments 7, 8, 10) and somewhat better 
for three -baseline designs with the exception of the extremely weak recovery in 
Experiment 4. Many -station designs (Experiments 5, 6, 9) improve the recovery 
especially in the equatorial configuration of Experiment 6. 

(e) £1 and ©o are the only parameters favored by the artificial equatorial 
design of Experiment 6 vs. the meridian one of Experiment 9. Recovery of the two 
components £ and 7 } of polar motion is not equal but depends on the position of the 
stations in each experiment with the exception of Experiment 6 . On the contrary, 
both the E and II parameters of precession-nutation are recovered with practically 
the same precision for every design. 
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The recovery of earth rotation parameters varies with the accuracy (Oq) 
of radio source coordinates . However, this dependence is much stronger for 
S, H; milder for the rest of the parameters; and even disappearing for £2, £, 77 in 
the strongest designs of Experiments 6 and 9 . 

(f) The results of the idealized designs of Experiments 6 and 9 can be used 
to verify some of the results of Chapter 6. Recovery of £, 77 , E, H is better for 
Experiment 9 than 6, while recovery of © 0 , Cl is better for Experiment 6. This is 
explained by the fact that Experiment 9 involves baselines parallel to the equator 
(16-17, 19-20) observing polar radio sources, a design sensitive to £, 77 , E, H. 
Experiment 6 involves equatorial baselines observing equatorial radio sources, the 
design most sensitive to 0 O , Cl (see Table 6. 1). 

(g) The recovery of baselines varies with station configuration. However, 
larger number of baselines does not in general improve the individual recovery as 
in the case of earth rotation parameters. For example, recovery in the three- 
baseline design of Experiment 2 is much better than in any of the many-station designs 
(Experiments 5, 6, 9). The relative recovery of angles in the same experiment is 
inversely proportional to the recovery of the opposite baselines. Dependence on 
radio source accuracy varies and disappears in the case of the equatorial configura- 
tion of Experiment 6, especially in comparison to the meridian configuration of 
Experiment 9. 

(2) Earth Rotation and Baseline Parameters As a Function of the Length 
of Observations and Earth Step 

For’ the variation of the recovery of earth rotation and baseline parameters 
with respect to the total time interval of observations T, it is noted that the standard 
deviations cr of the parameters are directly proportional to the standard deviations Ob 
of the observations. (Note that cr q = cr b . ) For this reason only the ratio p = a/oi, is 
depicted in Figs .9.4 and 9.5. 

The dependence of p on the earth step At is very strong for earth rotation 
parameters, especially for Cl (see Fig. 9.4), while not significant for baseline 
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parameters (see Fig. 9.5). In general, the recovery improves as At increases 
since the corresponding number of unknown parameters in the adjustment decreases. 
For earth rotation parameters (Fig. 9.4), standard deviation decreases strongly as 
T increases from 1 to 5 days. Anadditional increase of T from 5 to 10 days offers 
a small improvement, while no significant additional improvement is gained from 
the extension of T from 10 to SO days. The same is true for baseline parameters 
(lengths and angles, see Fig. 9.5) with the exception of an additional small improve- 
ment with the increase of T from 10 to 30 days. 

(3) Earth Rotation and Baseline Parameters As a Function of Earth Step 

As mentioned, the recovery of earth rotation and network geometry parameters 
improves with larger earth step At (Figs. 9. 6 and 9. 7). The only exceptions are with 
respect to the £,77 parameters in Experiments 1, 2, 4, 10 and with respect to 0 O in 
Experiment 4 where At = 24 hours results in worse recovery than At = 12 hours or 
even At - 6 hours in some cases. This "inversion" is correlated to the weak recovery 
of polar motion in Experiments 1, 2, 4, 10 (see Fig. 9.2) and to the weak recovery of 
©o in Experiment 4. On the contrary, for Experiments 6 and 9, where £, r\ and ©o 
recovery is better, variation of p with respect to At is "normal. " Such variation is 
negligible for © 0 in Experiment 6 and also for network geometry parameters in 
Experiments 6 and 9. Naturally an increase in the length of the earth step is a dis- 
advantage when short periodic variations in the earth rotation parameters are sought. 
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Fig. 9.2 Recovery of earth rotation parameters, Experiments 1 - 10. 


Total time span of observations: T = 1 day 
Earth step: At = 1 day 

a = a posteriori standard deviation of earth rotation parameters 
Oq = a priori standard deviation of radio source coordinates 
CTb = standard deviation of observational errors 

Note: Angular quantities when expressed in length units are 
arc lengths on a sphere of radius R = 6 371 000 m 
(1" « 30 m). 
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Fig. 9.3 Recovery of baseline lengths and angles, Experiments 1-10 


Total time span of observations: T = 1 day 
Earth step: At = 1 day 

For explanation of other symbols, see Fig. 9.2, p. 40. 
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Fig. 9 . 3 Recovery of baseline lengths and angles . 
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EXPERIMENT 5 

Fig. 9.3 (cont'd) 
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Fig. 9.4 Variation of earth rotation parameter recovery with 
stepsize and the total time interval of observations 
(Experiment 2). 
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stepsize and the total time interval of observations 
(Experiment 2). 
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Fig. 9 .6 Variation of earth rotation parameter recovery 
' with earth step. 
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APPENDIX A 


Inner Constraints and Their Relation to the "Geometry" of the 
Pseudoinverse of an Operator 

Suppose that a system of linear equations 

N x = y 

uXu u x 1 uXl 

is given (N, y known; x unknown) which is consistent (i.e. , it has at least one solution), 
and furthermore it has more than one solution. The objective is to show that there 
exists a set of linear constraints 

E T x = 0 

called inner constraints, such that the solution to the simultaneous system of 
equations N x = y and E x = 0 is unique and identical to the solution Xo = N y, 
where N + is the pseudoinverse of the matrix N. 

The matrix N is a representation of an operator N: X Y, where X and Y 
are both u-dimensional inner product spaces with elements u X 1 vectors and the 
usual Euclidean inner product 

< f, g> = g T f f, gEX or f, gEY 

Two subspaces can be introduced with respect to the operator N, the range R(N) = R 
of N defined as 

R - [y; y£Y and y = Nx for some xEX], RC Y 
and the kernel or null space K(N) = KCX of N, 

K = [x; xEX and Nx = 0] 

The consistency of the linear equations Nx = y refers to the fact that yGR, while 
the non-umqueness of solution to the fact that K Pi {o} c is non-empty, i.e. , K has 
other elements in addition to 0. If x' is a solution to Nx = y, then x =x' + x is also 
a solution where x is any element of K. Wc can therefore identify the solution 
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space S y C X, defined as 


S y = [x; xGX and Nx = y, y E R fixed] 

with a linear variety 

S y = [x; x = x' + x, x E K, x' fixed with Nx' = y] 

If K x is the orthogonal complement of K with respect to X, the projection xo = ^-1 (x) 
of any element x E Sy on K x is unique, and furthermore [Dermanis, 1977, §3.5] ‘ 

|| Xo|| = min || x |] 

XESy 

We can now define the pseudo inverse N + of th e operator N as .an operator 
N + : R=> X where 

N + y = x 0 = ^x (x) for any y E R and x E S y 

N + is therefore the ordinary inverse of the restriction of the operator N to K x (see 
[Desoer and Whalen, 1963, p. 444] for a more rigorous definition). The domain of 

-j- 4- 

definition of N can be extended to the whole space Y (N : Y X) by setting 
N + (y) = N + (y') 

where y‘ - (y) is the projection of yEY on R. 

The range of N + is K x , while its kernel (null) space is R x . For yE R we have 
Xo = N + y and Xo can be alternatively defined as the unique element in the intersection 
Sy D K x . This means that Xo is uniquely defined by its two properties: 

Xq E S y and xp E K x 

The first property is expressed by the original equations Nx = y, while a similar 
expression must be established for the second property. Let r = rank (N) < u and 
s = u - r be the rank deficiency of N . The dimension of the linear subspace K C X 
is then s, and suppose that tei], i = 1, 2, . . s is a basis in K. Since xp E K X , 
we have XqI ei for i - 1, 2, . . . , s, i e. , 

< Xo, ei > = el xp - 0 for 1=1, 2, . . . , s 
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Setting 


E = [% % ... e 6 ] 
u X s 

the above equations can be written in the compact matrix form 
E t Xo = 0 

Therefore, xo is uniquely determined as the vector satisfying both equations 
Nx = y and E T x = 0 

where the columns of E constitute a basis in K, the kernel (null space) of N. The 
constraints E T x = 0 are referred to as the "inner” constraints. 

The problem has now been reduced to that of finding a basis in K, i.e. , of 
finding a set of s linearly independent u X 1 vectors [ei] satisfying 

N et = 0 fori=l, 2, ..., s 

The equations above can be written again in matrix form 

N E = 0 where E = [ei eg ... e s ] 

We therefore only need to find a u x s matrix E with rank (E) = s = u - rank (N) 

(to secure linear independence of its columns ei) satisfying NE = 0.— Summarizing, 
we have the following theorem 

Given the system of linear equations 

N x = y (N, y known; x unknown) 

uXuXl uXl 

where rank (N) = r<u, s = u - r , and a u X s matrix E 
with rank (E) = s satisfying N E = 0, the vector x satis- 
fying simultaneously Nx = y and E T x = 0 is unique and 
identical to x = N + y, where N + is the pseudo inverse of N. 
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APPENDIX B 


Algorithm for First-Order Partitioned. Linear Regression 
Including Inner Constraints 


The problem in question is the construction of an algorithm for the least 
squares adjustment of a group of linearized observation equations of the form 

Vi = AiXi + AiX + Lt i = l, 2, m (A 1) 

with E CVi} = 0 and E{Vi vll = Pi. X refers to parameters common to all sets 
of observations, while X t refers to those appearing only in the i ul set of observations. 
Each set of observation equations gives rise to a contribution to the normal equations 
of the form 

Ni Ni 
Ni T Ni 

where 

N t = A! Pi Ai Ut = Ai Pi Li 

Ni = A* Pi Ai Ui = Aj T Pi Li 

and 

Ni = A/ Pt Ai 



Xi 


Ui 



+ 



X 


Ui 


(A 2) 


The combination of such individual sets of normal equations leads to the final set of 
normal equations incorporating all observations 


’n it" 


X 


u 




+ 


__ N t . N. 


.X 


L u_ 


where 



" 



* 

- 

Nx 

0 

0 


Nx 


Ux 

0 

Ns 

, N = 

Ns 

, U = 

Us 

0 

n b 



- 

U D 


X = 


Xi 

x 2 

x a 
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N 


= £ Nj , U = £ Uj 

i 1 


In case N = 

be solved through a First-Order Partitioned Linear Regression scheme as explained 
in [Brown and Trotter, 1969; and Uotila, 1973], We give here an extension of this 
scheme for the case when N is singular and a unique solution is ascertained only 
after the introduction of a set of inner constraints 

E t X + E t X = 0 (A 4) 


N N 
N 


is nonsingular (det N r 0), the normal equations above can 


with 

E t = [Ej ES ... El] = [Eo El 


E l] 


The normal equations can now be augmented as follows [Uotila, 1967]: 


- 


" tt - 


“ ^ - 

N N E 


X 


u 

T 




• 

N N E 


X 

+ 

u 

E T E t 0 


_ k l_ 


0 


(A 5) 


The new augmented coefficient matrix is nonsingular and inversion with the 
introduction of some new notation leads to 


"n 

N 

e" 

-1 


“ 

-1 

.. 




G 

G 

f" 





M 

M 


Q 

Q 





n t 

N 

E 






\ 

- 


G 

F 

_e t 

E t 

0 . 


M t 

M 


q t 

Q 


_p T 

f t 

F_ 


With the help of this inverse, the solution of the normal equations becomes 


-X 

= G U + 

G U 

-X 

= g t u + 

G U 


(A 7) 


If Gi, Gij denote submatrices of G and G, respectively, with dimensions equal to 
those of Ni and Nij, respectively, we have 


-X = £ G] Uj + G U 

(A 8) 

-X: = £ G:j Uj + G: U 
J 
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The inversion of Eq. (A 6) leads to 


Q = 

[M - 

M t M 1 Mf 1 


Q = 

- IT 1 

M Q 

(A 9) 

Q = 

M" 1 + 

M' 1 M Q M t M' 1 



Replacing M, M, M from (A 6) in terms of submatrices of N, N, N, and introducing 
the following notation 


Hi = N? N t 

Ri 

- Ni Ht 

Si = Nt - Ri 

S 

— £ Si 

H = £ Hi 

i 

K 

- £ N? 

i 

Li = Hi Ut 

L 

= £ Li 

Ti = Ni 1 Ui 

T 

= £ Ti 


we finally obtain after some algebraic manipulations 

S (E - H t Eo) 

Q' 1 = (A 10) 

_(E t - EjH) (-Ej K Eo) 

X = G L + F Eo T - G U (All) 

Xi = -T t ~ Ni [Eo F T (-U + L) + (Eo F Ej) T] ~ Hj X (A 12) 

G*j = 8i } Ni 1 + Hi G Hj + Hi (F Ej) N* + 

+ N? (Eo F T ) Hj + N't (E 0 F Eo) N* (A 13) 
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The algorithm can now be summarized in the following steps: 


1. 

Compute 

Ai, Aj, 

Li; E, 

i Eo 



2. 

Compute 

Ni, Ni, 

Ni, Ui, 

Ui 



3. 

Compute 

Hi, Ri, 

Si, Ti, 

Li 



4. 

Compute 

H, S, N 

, U, K, 

L, T 



5. 

Obtain G 

, F, F 

by inversion 




S 

E - 

H T Eo 

-i 

G 

F 

(E 

- H T Eo) T 

- Eo 

1 

o 


f t 

F 


6. Compute X from (All), Xi from (A 12), and Gij from (A 13) 
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APPENDIX C 


Computer Programs 

The simulation of observations and their adjustment are carried out 
with the help of three programs, VLBI SIMULATOR, PROG01, and PBOG02 
as explained in Chapter 7. These programs are presented here in detail together 
with their supporting subroutines . 

1 . Name of Program: VLBI SIMULATOR 

Function : Simulates distance and distance rate VLBI observations corres- 
ponding to time delays and their derivatives, according to a previously 
decided observational pattern. 

Input : Number of participating stations and radio sources, radio source 
and station coordinates, initial and final epochs of observations, time interval 
between successive observations, identifying numbers of stations and radio 
source participating in each observation. 

Output : Identification numbers for participating stations and radio source, 
observed (perfect) distance and distance rate, epoch of observation. 

Subroutines required: 

GRESID: Calculates Greenwich sidereal time for a given epoch 
JULIA: Converts Universal Time to Julian date 

EARTH: Provides simulated earth rotation parameters for a given epoch 
REDPI: Converts an angle in radians to the (- tt , ff ) interval. 

The following is a listing of VLBI SIMULATOR program and supporting 
subroutines. 
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RELEASE 2.0 


MAIN 


DATE = 76300 


21/57/44 


0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 
OOlo 

0u20 


0021 

CG22 

0023 

0024 
0 02 5 

0026 

CC27 

0028 

002 ° 

OwSO 

0031 

0032 
003 3 

0034 

0035 


003 6 


0037 

0038 
Ow39 


0 040 


0041 

0042 

0043 

0044 


0045 
0 046 


O****** CONVERT RA AND DEC TO RADIANS 
RAR=RA(I)*PI/1PO.DO 
DECR=D£C CI)*PI/1S0.D0 
uP=DCOS(DECR) 

SD=DS IN(PECR) 

CA=ncOMRAO) 

SA=DSIN(RAR) 

F1(I ) =CP*C A 
E? II )=CD*SA 
F3 < T )-SD 
5 CONTINUE 
C 

c 

c ******* READ STATION COORDINATES 

C*4****v PHI (LATITUDE) AND LONG ARE IN DEGREES, R IS IN METERS 
DO 6 1=1, IN 

Rf AD (5 ,701 ) IDUM,IL1,1L2,L0NG(I),PHI(I) ,R(I) 

7(>1 FORMAT ( 1X,I2,1X,A4,A2,3F15.1) 

WRITr 16,701 ) I,ILl,IL2,LDNG(I),PHim,R(I) 
loi FORMAT (3FH. 5 ) 

C*«***A CONVERT PHI AND LONG TO RADIANS 
PHIR=PHI (I )*PI/180.D0 
LONGR=LONG E I ) *PI/1 80 . DO 
C F =DCOS ( PH IP I 
SF=nsiN(PHIR) 

cl=dcos eloncr ) QP 

SL=ORIN(LQNGR) “g 

y(D = CF*CL 
Y 1 1 ) = CF $SL 
Z 1 1 ) = S F 

A CONTINUE 
C 
C 

Of G= p I / 12 . DO 
C 
C 

C******* R c AD INITT AL EPOCH IN UT K CO 

READ (5,201) IYEAR, IMO, id/i y.IHDUR.IMIN, SEC 
Wf\TT c (6,?Pl ) IYEAR,IMO,IE>AY, I HOUR ,IMIN,S C C 
201 FORMAT (1X,14,4T5,F15.1) 

C 

C 

C ******* cjno GRFEt 'W ITCH SIDERIAL TIM* AT EPOCH TO 

CALL GRESIDE IYEAR, IMO , IP A Y , IHPUR, IMI N, SEC.THO) 

C 

C 

C******* pf A q pinal EPOCH TF AND TIMF INTERVAL INHOURS 
READ! 5,201) JY fcAR, JMO , JO AY, JHOUR, JMI N, S EC J 
WFITF (6,201) JY6AR,JMn,J0AY, JHOUR, JMIN, SEC J 
RE AD ( 5 , 101 ) DT 
WK 1TC < 6, 101 ) DT 
C 
C 

IF(JYEAR.NF.IYFAR) GO TO 88 
IFf Jf.D.NF.IMO) GO TO 88 

C ******* NOTItr THAT IF TO, TF APE CLOSE RUT IN otFFFRENT MONTHS SECOND 
L ******* d/, t T TF MUST SE FXPRE C SFD AS A DAY O c PREVIOUS MONTH, f.G. AUG 33 
70 TF=D c LOAT(jr‘AY-IE'AY) > !'?4,riO+OFLOAT(J lJ nUR-IHOUR)+OFLOAT(JMIN-IMIN) 
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RELEASE 2.0 


HAIM 


DATE = 76300 


21/57/44 


0048 

0049 

0050 

0051 

0052 

0053 

0054 
0065 

0056 

0057 

005C 

0059 

0060 
0061 
006? 

0063 

0064 

0065 

0066 
C067 
0068 

0069 

0070 


7/60. D0+! SECJ—SEC )/36GO ,D0 
C 
C 

T0=0-D0 

T=0.00 

33 IF(T.GT.TF) GO TO 40 
C 

DO a 1=1 » IH 

c*****4* TRANSFORM QUASAR UNIT VECTOR TO EARTH FIXED SYSTEM AT EPOCH T 
ANCLP=OMG*(T-TO)+THO 
C A=DCOS (ANGLE) 

SA=DSIN(AMGLE) 

01=CA*E1(I)+SA*E2!I) 

Q2=-SA*F1{I)+CA*E2!I) 

Q3=E3!I) 

C 

DO a J = l» IN 

ZNT=PARCOS IX ( J )*Q1+Y ( J >*Q2+Z( J)*Q3) 

A RGU=X < U ) *Q 1 +Y I J ) *Q2+Z I J ) *Q3 
IF(ARGU.LT.-l.DO) ARGU=-1. DO 
IF(ARGU.GT.l.DO) ARGU=1.D0 
ZNT=DARCOS ! ARGU) 

ZNT=DAPS(ZNT) *180. DO/PI 
ZNTMAX*60.D0 

IF (ZNT .GT. ZNTMAX) ZNT*10 .D0**10 
INrJL'X(I,J)=ZNT 
8 CONTINUE 
C 

WRITE I 6 »207 ) T 


0071 


207 


FORMAT ( * 1 *//10X* ’EPOCH 


’ tF 15.5//20X»'1*»2X» '2 * ,2X, *3* t 2X, 


0072 

0073 


?*4',?Xt *5' ,2X, *6',2X, »7» t 2X,*8« r 2X, *9 * , 1 X, • 10* , IX , » 11 • f IX, *12 • , 
?1X, • 13*, IX, • 14 », IX, *15 SIX, *16' ,1X» »17S1X»*18',1X,'19»,1X,*20», 


?1X,'21*,1X,*22’,1X,*23*,1X,*24*,1X,»25*,1X,»26S1X, *27 SIX, *28*, 
?1X, *29 *, IX , *30 •//) 

C 


0074 


DO 9 1=1, IM 

0075 


WRITE (6. 10 5) I,RA(I),DEC!I ) , ( INDEX! I , J) , J=1 ,IN ) 

0076 

105 

FORMAT! IX, 12, 2F7.0, IX, 30 I IX. 12)) 

0077 

9 

C 

CONTINUE 

0078 

C 

INM=IN-1 

0079 


DO 555 1=1, IM 

0080 


00 555’ J=1 , I NM 

0081 


I F ( INDEX! I ,J ) .GT .1000) GO TO 555 

0C82 


JP=J4l 

ooe3 


DO 555 K=JP, IN ' 



0084 


IFITNOEXt I,K J.GT.1000) GO TO 555 

OOP 5 


WRITEC7 t ?05) JtK,I t T 

0086 

705 

FORMAT ( 3 1 5 t PI 5 .2) 

0087 

555 

C 

CONTINUE 

0088 

c 

99 

T=T+DT 

00 89 


GO TO 33 

0090 

88 

WHITER *500) 

C0°1 

500 

FORMAT ( ' 1 ' /// ////lOX* • WRONG INITIAL AND FINAL DATE'/ 
?10X» • PROGRAM INTERRUPTED THROUGH STATEMENT 88») 

0092 

40 

STOP 

009 3 


END 
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0008 

0009 

0010 

oou 

0012 

0013 

0014 

0015 

0016 

0017 

0018 
0C1O 
0020 
0021 


0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 
0022 

0033 

0034 

0035 


0026 

0037 

0038 

0039 


0040 

0041 
004? 


0043 

0044 

0045 

0046 

0047 
00^8 


C******* ra AND DEC ARE IN DFGREES 
DO 5 1=1, IM 

RE AD (5, 321 ) IDM,IL1»IL2,IL3»RA(I) , DEC (II 
WRITE (6,321 1 I,IL1,IL?,IL3,RA(I) ,DEC(I) 

321 FORMAT (lXtI3,7X»3A4»F17.7»F2Q.9) 

C ******* CONVERT RA AND DEC TO RADIANS 
RAR=RA(I)*PI/1PO.DO 
DECR=D6C ( I )*PI/180.D0 
CD=DCOS (DECR > 

SD=DS I N (DECR ) 

CA=DCOS (RAR) 

SA=DSIN(RAR) 

Eld )=CD*CA 
E 2 ( I ) =CC'4S A 
E3 ( 1 1 =SD 

5 CONTINUE 
C 

c 

C ******* READ station coordinates 

c******* phklatitudei and long are in degrees, r is in meters 

DO 6 1=1, IN 

READ (5, 701) IDM,ILltIL2,L0NG(I),PHI (I ),R(I) 

WRITE (6, 701) I,IL1,IL2,L0NG(I),PHI(I),R(I) 

701 F0RMATUX,I2,1X,A4,A2,3F15.1) 

C******* CONVERT PHI AND LONG TO RADIANS 
PHIR=PHI ( I)4PI/180.D0 
LONGK=LONG ( I J *PI/18Q. DO 
tF=DCOS(PHIR) 

SF=DS IN ( PHIR ) 

CL=DCOS (LGNGR ) 

SL=DSIN( LONGR) 

X( I )=CF*CL 
Y ( I ) =CF*SL 
zr I)=SF 

6 CONTINUE 
C 

C 

C******* READ CORRECTION TO GREENWITCH SIDERIAL TIME DDTH IN ARC SEC 
READ (5,101) DDTH 
WR IT S (6 , 1 01 ) DDTH 
101 FORMAT ( 3F15 .5 ) 

C******* CONVERT TO RADIANS 

DDTH=DOTH*PI/( 180.00*3600.00) 

C 

C 

C ******* READ INITIAL EPOCH IN UT 

RE AD (5, 201) IYEAR,IMO,IDAY,IHOUR,IMIN,SEC 
WRITE <6, 201 I IYEAR,IMO,IDAY,IHOUR,IMIN,SEC 
201 FORMAT (5I5,F15.5) 

C 

C 

C******* PINO GREENWITCH SIDERIAL TIME AT EPOCH TO 

CALL GRESIDflYEAR , I MO, ID AY ,1 HOUR » IMIN ,SSC,THO ) 

COPY( 1)=TH0 
THO=TH04(lpn. DO/PI) 

WRITE ( 6 , 136 ) THO 

136 FORMAT(///10XV'APPR. THO= SF25.9//) 

WRITE (7 , 101 ) THO 
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0002 
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0007 


RELEASE 2.0 


MA IN 


DATE = 76300 2 1/57/12 


C*****************************************************:*********:i ****** ********** 
C********* 1 ********** 

C ********* VLB1 NETMORK ♦ OBSERVATION SIMULATOR f ********** 

Q********* | ********** 

C * ** ******* ** ******* **** **** ******* ******** **************************** ********* 







CONTINUOUS 

VERSION : EARTH MOTION PARAMETERS ********** 



PROVIDED BY SUBROUTINE EARTH 







c ***** ***** 



******❖->** 

c********* 

THE MODEL 

USED IS: 

A#**####** 

£********* 




£********* 


T j 

*#*##****# 

C********* 

0(1, J,P,K) 

=<XJ-XI) *S(XP,HP)* 

*#*#*#**** 

£;*****> **** 


T 

**$******* 

C********* 


*R3 { OMG*TK +THO ) *S <XX,HH)*QP 

*.#$******* 

c ->** ****** 


1 

********** 

C********* 

( 

1 0 A ) 

********** 

C ****** * A* 

S( A, 3) = { 

0 1 -B ) 

********** 

£********* 

( 

-A B 1 ) 

********** 

C********* 



******* * ** 

C********* 

T 


********** 

£********* 

QP = (COS(DECP)*COStRAP) , 

♦ * * ** * * * * * 

C#** ****** 

COS(DECP)*SIN(RAP), SIN(DECP) ) 

***** ***** 

C********* 



********** 

C****M**ft 

XJ,XI 

: POSITION VECTORS OF J AND I 

*•* ******** 

£********* 


STATIONS IN FARTH FIXED SYSTEM 

********** 

c********* 

DECP ,RAP 

: DECLINATION AND RIGHT ASCENSION 

* ********* 

O******** 


OF P QUASAR IN Q'JASOR FIXED SYSTEM ********** 

£********* 

QP 

: UNIT VECTOR OF P QUASAR IN 

********** 

C********* 


QUASAR FIXED SYSTEM 

********** 

£********* 

OMG 

: ROTATIONAL VELOCITY OF EARTH 

********** 

Q********* 

XP, HP 

: POLAR MOTION SMALL ANGLES 

** ***** 

C********* 

XX, HH 

: PRECESSION NUTATION SMALL ANGLES 

********** 

C********* 

D(I,J,P,K) 

: OBSERVED DISTANCE INVOLVING IJ 

********** 

C********* 


BASEL INF AND P QUASAR AT TIMF TK 

********** 

Q********* 


AFTER INITIAL REFERENCE EPOCH TO 

********** 

£********* 

THO 

: GRFENWITCH SIDERIAL TIMF AT FPOCH TO ********** 

C********* 

THFTA 

t UT1 ANGLE 

* * * * * * * s\ * lie 

C********* 


r 

* ** * * * s{r * * 



C **** ***** ********** 

IMPLICIT REAL*8(A-H,L-Z) 

DIMENSION X(4) ,Y(4) ,Z(4),El(l8> ,E2(18), E3 C 1 8 ) , RA ( 1 6 ) ,D EC ( 18) 

2 » PHI (4) ,L0NG(4),R(4),C0PY(3) 

C X(INJ,Y(IN) ,Z< IN )»E1(IM),E2(IM),E3(IM), INDEX! IM , IN) , RA ( IM ) ,D£C (IM > , 

C PHI ( IN) ,L0NG( IN ) fR(IN) 

C I N=N0 OF STATIONS, IM=N0 OF QUASARS, IN.LT.35 

C 

PI=4.C0*DATAN( 1.D0) 

RC ART =6371000. DO 
C 

READ { 5 , 100 ) IN, IM 
WRITE ( 6 , 100 ) I N, IM 
100 FOkMA T ( 21 5 ) 

C 

C ******* READ QUASAR COORDINATES 
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MAIN 
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0049 


0050 

0051 

0052 
0063 


0064 

0055 


0056 


0057 

0058 

0059 

0060 
0061 


0062 


0063 

0064 

0065 

0066 

0067 

0068 
0069 
0C70 
0071 
007? 

0073 

00 74 

0075 

0076 

0077 
OC78 

0079 

0080 

ooei 

*0082 

0083 

0084 

0085 

0086 
0087 


THO=COPY(1)+ODTH 

C 

c 

c**«***# READ final epoch tf and time interval inhours 

READ (5, 201) JYEAR, JM(1» JDAY ,JHOUR* JMIN.SECJ 
WRITE ( 6,201 I JYEAR, JKO, JOAY»JHQUR, JKIN,SECJ 
READ (5 f 101 ) DT 
WRITE ( 6r 101 ) DT 
C 
C 

IFUY6AR.NE.IY EAR) 60 TO 88 
IF(JMO.NE.IMO) GO TO 88 

C ******* NOTICE THAT IF TO, TF ARE CLOSE BUT IN DIFFERENT MONTHS SECOND 
C*****?* D aT 6 TF MUST BE EXPRESSED AS A DAY OF PREVIOUS MONTH, E.G. AUG 33 
79 TF -DFLOAT ( JO AY-10 AY)*24.00*DFL0AT < JHOUR-IHOUR ) +DFLDAT l JMIN-IMIN) 

?/60.D0 + ( SECJ-S EC )/3600.DO 

c 

c 

TO=O.DO 

T=O.0O 

33 IF(T.GT.TF) GO TO 40 
C 

555 READ (5, 107) I»J,K 

JFd.EQ.O) GO TO 99 

C ******* 1=0 IS A CODE TO INDICATE END OF OBSERVATIONS FOR THIS EPOCH 
C 

CALL EARTH (T,XP, HP, <3MG, THETA »XX»HK1 
C 

C MODIFY EARTH PARAMETERS AND WRITE / PUNCH 

XPOsXP+REART 
HPD=MP*RFART 
XXD=XX*R EAHT 
HHD=HH4RP*RT 
0MGO=OMG4( 12.DQ/PI ) 

THtT=THETA 

THET =THET* ( 180.D0/PI1 
ITH1=THET 

TH2= (THET-DFLOATt ITH1 ) )*60 .DO 
ITH2=TH2 

TH3=(TH2-DFLr)AT( ITH2) ) *60.00 

WRITE (6, 801) T,XPD,HPD»XXD,HHD,ITH1,ITH2,TH3.0MGD 
8 01 FORMAT 1 3X»F7 .2 ,4F10«4 , IX,2 I4,1X ,F8.4,F12.9) 

C 

COPY ( 1 )= E11K) -XX*E3(K I 

COPY ( 2 )= + E2 UO+KH4E31K ) 

COPY (3 ) = XX*£1 1 K)-HH#E2(K)+ E?(K) 

C 

'ANGLE=THETA+THO 

CA=DCOS(ANGLE) 

SA=DSIN(ANSLE ( 

c 

01= CA*COPYll )+SA*COi>YC2> 

02=-S A*COPYt 1 ) +C A*COPY ( 2 ) 

Q3=C0PY{ 3) 

c 

QFl=-SA*C0PYm+CA*CQPY(2) 

QF2=-CA*C0Prf 1 1 -$A*C0PY1 2 ) 

0P3=-XP*QFUHP*QF2 
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C 


0 088 


COPY! 1 >= 01 +XP*Q3 

0089 


C3PY ( 2 ) = Q2-HP*03 

0090 

r 

COPY{3)=-XP*01+HP*Q?+ 03 

t 

00°1 

w 

rs=<x( j>*r(j>-xii men )*coPY<n + (Y(jm (j>-Ym*R(in*coPYt 2 ) 


r 

? -Mzum{j)-zm*K(in*coPYm 

009 2 

V 

PRNG=(X( J)*R( J)-X<n*R(I) )*QFl + (YtJ)*R( J)-Y{I)*Rm )*QF2 
2+(Z[J)#R<J)-Z(I)*R(I))*QF3 


C 

FRSG= DERIVATIVE OF DELAY DISTANCE OS ( METERS/HOUR ) 

0 093 

r 

FRNGsFRNL^OMG 

0094 

L 

WRITE C 6 * 107 J I,J,K,DS,.PRNG,T 

0C95 


WRITF <7* 107) I,J,K,DS,FRNG,T 

0096 

107 

F ORMAT (315 »3 FI 5. 5) 

0097 


CO TH 555 

0098 

99 

T =T+DT 

0099 


GO TCI 33 

C 100 

88 

WRITC(6,500) 

0101 

5 00 

FORMAT (’1*//// ///l OX* * WRONG INITIAL AND FINAL DATE*/ 
71 OX, 'PROGRAM INTERRUPTED THROUGH STATEMENT 88*) 

0102 

40 

STOP 

0103 


L-ND 
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GRESID 


DATE = 76300 


21/57/12 


0601 


OC 02 
0 603 
O-'OA 
GC05 
CC06 

0007 

0008 
nr09 
OCIO 

ocai 
0 012 
0613 
0 016 

0015 

0016 
0C1 7 
0018 


'U8R0UT1NE GRESID(IYFAR,IMO,IDAY,IHOUR,IMIN,SEC,THO) 

(,********* ^ #**#++#**%: 

£***************************, ********* **************************]£**************** 
£********* . ********** 

C.**** :**#** CALCULATES THE GREENWITCH SIOERIAL TIME ********** 

C ********* THO, FOR EPOCH WITH DATE 1YEAR IMO IDAY ********** 

C*****v*** IHOUR IMIN SEC UT. ********** 

C ********* ********** 

C ************************************************ ******************************* 

C********* ********** 


IMPLICIT REAL’S'Ct A-H,L-Z) 
FI=6.n0=t‘0ATAN( 1.00} 

11=0 
I 2=0 


A =0.06 


CALL JULIA (I YEAR, IMO, IDA Y , 1 1 , 12, A ,M JD ) 
T= (MJD- 26 1 6620 .DO) /36525 .DO 
T 2 =T *T 


TP0=6>V .690983210+3^000 .76S9DO*T +0 ,0003 8708DO#T2 

THO=THO*PI/lfcO.DO 

CALL REDPI (THO ) 

6T=DFLDAT( IHOUR) #60.D0+DFLOAT(IMIN)+SEC/60.D0 
PTHDT = 6 .375269 500/(10.00**3) 

THO=THO+PT*DTHDT 
CALL REOPI(THO) 

K-TJRN 

END 
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JULIA 


DATE = 76300 


21/57/12 


0001 


0002 

0003 

0004 

0005 

0006 
C007 
GOOF 
000® 
0010 
0011 
0012 

0013 

0014 

0015 


SUBROUTINE JULIAlIYFARflM, I0AY, IHHH , IMMM ,S ,MJD) ( 


£ *****:! *** 

i 

i 

******** ** 

C ********* 

1 

***> ii****** 

£********* 

SUBROUTINE TO CONVERT UNIVERSAL TIME TO i 

********** 

£******■! >** 

JULIAN DATE 

********** 


1 

********** 

c********* 

1 

A********** 


c **#*****$ A ****** *#:* ******* »**»*#*» &#j**;&*&>»****i*i)<**^:([* ******* *4*^ ****** ********** 
IMPLICIT REALMS {A-H*L-Z) 1 

DIMENSION I’:0NTH(12) ! 

DAIA IMQNTH /0»31,59t90*120*151»181 »212 »243 »273 »304» 334/ ! 

H=OFLOATt IHHH) j 

K=[iFLOAT(IMMM) 

ICHO=0 

101 5 = ( IYEAR-18®7)/4 ' 

I t «TYFAR.GT.l°00) IDT S=IDIS— 1 ' 

ICK=4*(IYEAR/4J , 

IFaYEAR.EO.lCH.ANP.IM.GT. 21 IC0D=1 
IF(IVEAR.£Q.1900J 1COD=0 

NJ0=241 5020+ H YEAR-1900) *365. D0+IDIS+IM0NTHUM)+IC0D-0.5d6+IDAY 
X+H/24 .D0+M/1440.D0+S/86400 .00 
RETURN 
l NO 
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000 ? 

0003 

0004 


oros 
0006 
C 007 

ooo a 

0009 

0010 
con 
0012 
0013 
00X4 

0015 

0016 
0 017 
00 If. 

0019 

0020 
0021 
0022 
0 023 

0024 

0025 

0026 
002 7 
002 8 
GC29 

0030 

0031 
0 032 

0033 

0034 
0C3 5 


IV 0 1 
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EARTH 


DATE = 76300 


21/57/12 


C 

C 

c 

c 

c 

c 


SUBROUTINE E ARTH IT ,XP ,HP ,OMG, THETA, XX, HH ) 

T : TIME INTERVAL AFTER INITIAL EPOCH (HOURS) 

XP,HP: POLAR MOTION COORDINATES (RADIANS) 

THETA: UT1 ANGLE IN RADIANS 

OMS : ROTATIONAL VELOCITY OF EARTH (RAOTANS PER HOUR) 

SMAUL ANGLES {RADIANS> ' 

DIMENSION P<10) ,A< 10) ,8( 10),PS(10),FC10) 

DATA P/8760. DO, 438 0.00,2190. DO, 720. DO, 360. DO, 

2165.00. 24.00.12.00.1 .DO, O.SDO/, 

?A/1. 1600,0. 63 DO, 0.2P.D0, 0.1 2D0, 0.0600, 

40. 0300 .0. 00200. 0.00100.0. 00000.0.00000/, 
5R/0.°500,0.70DO,0.??DO, 0.1300,0.0500, 

60. 01PO,n. 00400 ,0.00 100 ,O.OOODO,O.OOODO/, 

7PS/3. 15 DO, 5. 2 BOO, 1 .3?DO,0.15DOtO.CDO ,0.99DO,2*0D0» 

86. 1400,2.75 DO, 4, 1300/, 

^ / ^n^ D ° ,2 * P5D0 ’ 6 ' 17D0,A *^ 50C,0 * 05D0 » 5 ’ 18D0 ’ 3 * 75D0 » 1 *24D0i 
70.3500,1.6700/ T 


C 


5 

C 


91=4.00*047 AN ( 1 , 0 0) 

P 12=2 • DO’>P I 
RE ARTs 63 7 1000 . DO 
M10=6 .00 
M20=0 .00 
PC=429. 00*24. DO 
0=0005(212*7/90 
S=DSINtPI?*T/FC) 
M1=C*M10-S*M20 
M 2=S*M10+C*M20 ’ 

OO 5 r=i,io 

AkG=(PI2*T )/P< I) 

Ml =M1 +A ( 1 ) *DCOS ( ARG+P S ( I ) ) 
M2=M2+ n (I )*nCOS(ARG+ F (I ) ) 
CONTI'-'UE 


XP=-M2/R6ART 

HP=-M1/REART 


XX=0. 0600+0. 01 rococos (2. 34D0+T/24 . 00 ) -0 .003*DC0S(0 ,100+T/i 2.D0 ) 

xx=xx/ReART’ 02I:,0 * ^ ' COS^3 ' llDO+T/24 * DO) ~ 0 * 005DO,^ ‘ DCOS<0 * 07R0 ^ T/12 * 00, 


HH=HK/RSART 

TD=?P51 00+T/24 .00 

A RC-1 =P I 2*TO/?65 .2500-0 .49000 


ARU 2 = 2 . rjO*Pl2*TD/3 65. 2500-0. 862 D 0 

THETA=O.OO176DO*(T/24.D0)+25.1D0*DSIN<ARGl)-9.2OO*0SIN(ARG2) 

FACTOR = PI2/(!64OOO0O.OO 

THFTA=P 12* (T/24.D0 )+THETA*FACTOR 

^C)MG=0 .00175D0/ 24 .00+ ( 25 .1DO*DCOS ( ARG1 ) — 1 8»4DO*OCQS ( ARG2) ) ’ 
(PI 2/8766 *DG ) f 


0 P C,= P I 2 /2 4 . D 0* F A CT 0 R*OMG 


RETURN 

FNO 
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REDPI 


DATE = 76300 


0001 


0002 

0003 

0004 

0005 
0 006 
0007 
0000 
0 C 09 
0010 
0 C 11 
0012 

0013 

0014 


21^56/32 

i 


SUPRGUT 1 NE REDPI (A) 


C ****** ***** *** ******* a********************************************************* 
£********* 

C********* - 

C********* 

L ********* 

C *********************************************************** ********************* 
C********* ] ********** 


REDUCES AN ANGLE A IN RADIANS TO THE 
f -Pit +PI ) INTERVAL 


********** 

********** 

********** 

$C$t 


IMPLICIT REAL*e(A-H,L-Z) 
PI= 4 . 00 * 0 ATAN( l.DO) 

PI 2 =PI* 2.00 
C 0 PY=A/PI 2 
I A=COPY 

IF ( A ,GE . 0 . DO ) A=A-CFLOAT (IA >*PI 2 
IF(A.LT.O.DO) A=A+DFL 0 AT(IA)*PI 2 
5 IPIA.GT.-PI. AND. A. LF.P I ) GO TO 8 

IFtA.LT.O.DO) A=A+PI 2 
IF (A.GT.O.DO) A=A-PI 2 
GO TO 5 
8 RETURN 

END 



2. Name of Program: PROG01 
Function: 

Adjustment of a set of observations with total time-span equal to both 
earth and station step. Earth rotation parameters and station coordinates 
are treated as constants over the total time span of observations . A number 
of weighting options for radio source coordinates are available. 

Input: 

Number of weighting options for radio source coordinates, standard 
deviations of observations and radio source coordinates, approximate 
values of parameters, identification numbers of participating stations and 
radio source, distance and distance rate observations. 

Output: 

Standard deviations of adjusted parameters . 

Subroutines required: 

REDPI. 

(Following is a listing of PROG01. ) 
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MAIN 


DATS = 76300 


33/55/2? 


C, *** * ** ******** * **** ***************** ******** ************************* ********** 


[****»»** PRQG01 ********** 

£n******** ********** 

C********* VLSI ' ADJUST KENT PROGRAM ********** 

C********* ********** 

C********* THIS VERSION PFRPORMS AN APJUSTMFNT OF ********** 

C**»****** ONE CAY OF OBSERVATIONS TREATING ALU ********** 

C********* PARAMETERS AS CONSTANTS OVER THF ONE DAY ********** 

C********* INTERVAL. ********** 

c *„******* AN OPTION IS AVAILABLE FOR APPLYING A ********** 

L ********* NUMBER OE DIFFERENT WEIGHTS ON RADIO ********** 

C ********* SOURCE COORDINATES ********** 

C********* ********** 

C ********* OBSERVAT IONS • DELAY DISTANCES AND ********** 

C ********* THEIR DERIVATIVES ********** 

C********* ********** 

C********* UNKNOWNS: EARTH MOTION PARAMETERS ********** 

C********* STATION C OUASAR COORDINATES ********** 

C********* CLOCK OFFSETS AND DRIFTS ********** 

C ********* ********** 

t********* SYSTEMS DEFINITION: ********** 

C********* OUASAR FIXED: WEIGHTED OUASAR COORDINATES ********** 

O******** EARTH FIXFD: BY MEANS OF INNFR CONSTRAINTS ********** 

C********* ********** 

C********* SUBROUTINES REQUIRED : REDPI ********** 

C********* ********** 

C********* ********** 

C ******************************* ************************************************ 
C, ********* ********** 

o********* ********** 

C 

C A ROW OF DE5IGN MATRIX CONTAINING PARTTALS OF DELAY 

C OBSERVATIONS WITH RESPECT TO PARAMETFRS 

C B ROW OF DESIGN MATRIX CONTAINING PARTI ALS OF DELAY 

C DFRIVATIVE OBSERVATIONS WITH RESPECT TO PARAMETERS 

C CC COEFFICIENT MATRIX OF INNER CONSTRAINTS 

C N COEFFICIENT MATPIX OF NORMAL FOUATIONS 

c U CONSTANT VECTOR OF NORMAL EQUATIONS 

C XX VECTOR OF UNKNOWNS (CORRECTIONS TO PARAMETERS) 

C THE ORDER OF THE PARAMETERS IS: 

C XP, HP POLAR MOTION ANGLES KSI AND ETA 

C XN, HN PRECESSION-NUTATION ANGLES (CAPITAL) KSI C ETA 

C TH THETA ZERO ANGLE (CAST AT INITIAL EPOCH) 

C OMG ANGULAR VELOCITY OF EARTH ROTATION 

C XII ),Y(I) »Z(I) 

C COORDINATES OF STATION I ( 1=1 ,? ». . , IN) 

C RA(J), D(J) 

C RIGHT ASCENSION AND DECLINATION OF 

C RADIO SOURCE J CJ=1,2,...,IM) 

C CLOCK OFFSET AND DRIFT AT STATION I (1=1,2, IN ) 


IMPLICIT REAL*8(A-H,L-Z) 

DIMENSION A( 37 ) »N( 45 ,45) ,X (3) ,Y (3 ) , Z (3 ) ,RA ( B ) ,D ( 8 ) ,U (37) 
•2,XX(37),H (45) ,L2(45) ,CC(8,37),00(3) ,DD (3) , B(37 > ,NKEEP (45,45 ) 

3, SPOUA( 10 ) jPRATIO ( 10) ,F ( 13 ) 

DIMENSION A(KK ),N(K8,KB) ,X(N) ,Y <N ) , Z (N) , RA( M ) ,D ( N ) ,U(KK } ,XX ( KK ) 
2,L1(K8) ,L2(K8) , CC(8,KK) ,00(N) ,DD(N) 
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C 

3* B (KK ) ,NKEEP (K8*K8 ) *SDQUA( 10) *PRATIO( 10) 


C 

r 

WHERE: N=NO OF STATIONSt M=NO OF QUASARS* K=6+3*N+2*M, KK=K+2*N, K8=KK+8 


c 

REAP EXPERIMENT IDENTIFICATION NO 

0003 


R EAO ( 5 , 100 ) I EXP 

000 A 


WRITE ( 6*977 ) ICXP 

0005 

977 

FORMAT! '1 V///10X, 'EXPERIMENT NO =',I5/10X 


r 

2* ' *********#*#********»,///// ) 


c 

READ NUMBER OF WEIGHTING OPTIONS 

0006 


RE AO ( 5 * 100 J IOPT 

0007 


WRITE (6 *978 1 IOPT 

0003 

978 

r 

FORMATUOX, 'REQUESTED NO OF WEIGHTING OPTIONS =*,I5/////) 


c 

RF AD A PRIORI STANDARD DEVIATIONS: 


c 

SD8D : OF DELAY DISTANCE IN CM 


c 

SDBF : OF DELAY DISTANCE DERIVATIVE IN M/HOUR 


c 

f 

SOQUAJ OF QUASAR COORDINATES IN CM 

0009 

V 

RE AO ( 5 * 101 ) SDBD*SDBF 

0010 

101 

r 

FORMAT! 5F15. 5) 

0011 

w 

DO 1 8 1 1=1, IOPT 

0012 


READ! 5*101) SDQUA(I) 

0013 


PRATIO!I)=SDBD**2/S0QUA( I)**2 

0014 

181 

C 

CONTINUE 

0015 

c 

R=6371000.D0 

0016 


CL IGHT=3.D0*( 10.00**8 > 

0017 


CLIGHT=CLIGHT/R 

0018 


PI=4.DO*OATAN(1.DO) 

001« 


rioo=r*ioo.do 

0 020 


RHKDM<?=!12. DO/PI ) * (10.00**9 ) 

0021 

r 

RA£M3=(1B0. DO/FI) *2600000. P0 

0022 

U 

PF={SDBD/R100)/(SDBF/R) 

0023 


PF=PF**2 


c 

WEIGHTS : 


c 

c 

PF IS NOW IN (EARTH RADII/H0UR)**-2, PRATIO IN (EARTH RADII) **-2 


c 

READ IN= NO OF STATIONS S IM= NO OF QUASARS 

0024 


RE AD ( 5 1 1 00 ) IN, IM 

0025 

100 

FORMAT! 515 ) 

0026 


WRI T E ( 6*920 ) IN*IM 

0027 

9 20 

FORMAT (5X, 'NO OF STATIONS ='*I5/5X»'N0 OF QUASARS =»,I5) 

0028 


K=6+3*IN+2*IM 

0C29 


kk=k+?*in 

0030 

r 

KS=KK+? 


w 

c 

INITIALIZE POLAR MOTION (XP,HP) AND PRECESSION-NUTATION ( XN* HN) 


c 

PAKAMETFRS ( APPROXIMATE VALUES) TO ZERO 

0031 


XP=0. DO 

0C32 


HP=O.DO 

0033 


X N=0. DO 

0 C34 

r 

HN=0. DO 

0035 


DO 51 I=1*IN 
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0036 

0037 

0038 


00i«» 

0040 

0041 

0042 

0043 

0044 

0045 


0046 

0047 
0C48 

0049 

0050 

0051 

0052 
0C53 

0054 

0055 

0056 

0057 

0058 


OC.59 

0060 

0061 
0Cfc2 
0063 
0 064 

0065 

0066 
0067 


0068 

0069 

0070 

0071 


Dom=o.oo 

DO (I )=0.D0 
51 CONTINUE 
C 

C ***** ********** *********** a.********************#******************************* 

c read approximate values of parameters theta zero, omg, 

C Xtl), Y(I), Z(I), <I=1,IN), RA(J), DCJ), (J=1,IM> 

C ANGLES ARE IN DEGREFS » DISTANCES IN METERS AND OMG IN RADIANS PER HOUR 

C ************ *****= ************************************************************** 

c 

READ < 5 * 101 ) TH 

WRITE (6 tC60) TH 

860 FORMAT {//5X, ’THETA ZERO = SF20.9/) 

TM=TH*P 1/180. DO 
OHG=P 1/12 • DO 
C 

WRITE J6.450) 

450 FORMAT (////♦ 10X. ’STATION 6 RADIO SOURCES APPROXIMATE COORDINATES’ 

7.//3X, 'STATIONS', 10X, 'LONG* 

Z.lZX.'LAT’.llX.'R’./l 

C 

DO 5 1*1, IN 

READ ( 5 , 701 ) TDM»ILl,IL2,AL0NGtAPHI,AR 
WR1TK« 6,701) I.IL1.IL2, ALCNG,APHI, AR 
701 FCRMAT{1X,I2»1X,A4,A2,3F15«1) 

APHI=APHI*PI/1P0.D0 
AL CiNG* ALON G* P I /I BO . DO 
C CONVfRT TO CARTESIAN 

XII }=AR*C'COS ( APHI 1 *DCOS ( ALONG 1 
Ym*AR*r , COS(APHI)* , OSIN(ALONG) 

Zm=AR*OSINCAPHI ) 

C CONV c RT TO EARTH RADII UNITS 

X( 1> = X( 1 l/R 
Y(I)=Y<I)/R 
Z(’)=Z( I)/R 

5 CONTINUE 
C 
C 

WRITE {6,451 1 

451 FORMAT (//l IX, ’RADIO SOURCE ', 10X, ’ RA 17X , ’D EC* ,//) 

C 

DO 6 J=1,IM 

REAO (5,321) IDM,IL1»IL2,1L3,RA(J),D(J) 

WRITE (6,321) J,1L1,1L2,IL3,RA(J),D(J) 

321 F0RMATUX,I3,7X,3A4,F17.7,F20.9) 

C CONVERT TO RADI AN UNITS 

RAU)=«A(JJ*PI/180.DO 
D(J) = D( J)*PI/ltO.DO 

6 CONTINUE 
C 

£**********+*******•************************************************************* 
C INITIATE N, U, XX, CC TO ZERO 

£ 3 $ 4:4c ******* ****** *********************************************************** 

C 

DO 7 1=1, KK 
XX { I ) =0 .DO 

um=o.oo 

DO 7 J=1 ,8 
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0072 

0073 

0074 

0075 

0076 

0077 
00 7 E 


0079 

0080 
0061 

0082 
008 3 

0084 

0085 

0086 
ooe? 
008 £ 
0089 
C090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 


0100 

0101 

0102 


0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 


CC(J,I)=0.D0 
7 CONTINUE 

C 

DO 16 1=1, K8 

00 16 J=l,K8 
N( I,J)=O.DO 
NKEEFM I »J)=0.D0 

16 CONTINUE 

C 

C FORMULATE CONSTRAINT MATRIX CC 

C 

CC (5, 1 ) =1 .DO 
CC(4,2)=I .DO 
CC ( 6, 5 ) =-l .DO 
C 

DO 3 1=1, IN 

1 61=9 + 3*1 
IB 2=1 8 1 + 1 
183=182+1 
CC(1, IB1 ) = 1 .DO 
>CC ( 2, 1 82 ) = 1 .DO 
CC (3,183 1 = 1. DO 
CC(4,IP2 )=-2tI) 

CC (4, I P3 J =Y ( I ) 

CC(5,iei)=2<I) 

CC (5, 1 83 ) = -X (I ) 

CC(6,IB1)=-Y(I) 

CC (6,182 )=X{ I ) 
ici=k+i 

IC2=IC 1+IN 
CC (7 , 1C 1 ) = 1 .DO 
CC(£,IC2)=1.D0 
3 CONTINUE 
C 

C****************#***************************************-**************** *****## 

c 

c READ OBSERVATIONS AND IDENTIFIERS: 

C DS= DISTANCE IN METERS OF IJ BASELINE OBSERVATION TO QUASAR IP 

C AT TIME TK (IN HOURS) AFTER SOME INITIAL EPOCH TO 

c 

C *****$*#$*$* *4 4 444 * ./rfr*#********^******##*^* ******!>*******)<******:***** ********** 

C 

c 

IC0UNT=0 

^ i 

33 . READ ( 5 , 104 ) I , U, IP ,DS , FRNS »TK , ICHEK 
104 FORMAT! 3 15 .3815.5,215 ) 

C 

C ICHEK IS A CODE INDICATING END OF OATA FOR ICHEK=1 

IF(ICHEK.EQ.l) GO TO 66 
IC0U\T=IC0UNT+1 
WRITE ( P , 100) I,J,IP 
C 

96 DS=DS/R 

FRNG=FRNG/R 
DX=X( J ) -X { I) 

D Y=Y { J ) —V C I > 

DZ=Z(J)-Z(I) 



MAIN 


DATS = 76300 


C'S*- 

FORTRAN IV G1 RELEASE 2.0 


0111 

0112 

0113 

0114 

0115 

0116 
0117 
CUE 

0119 

0120 


0121 

0122 

0122 

0124 

0125 

0126 
0127 


0128 


0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 
0 145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 


CA=DCOS ( RA ( IP ) ) 

$A=n$IN(RA(IP) ) 

CP=DCOS (D( IP 1 ) 

$0=05 1 N (0 ( TP ) ) 

Y1 =OMG*TK+TM 
Y2=Y1-RA{ IF) 

CK=DCC<;fYll 
SK=0 S IN l Y1 ) 

CKF=rC0S(Y2) 

SKP=DSIN t Y? ) 

C 

C FORMULATE PARTIALS 

C 

Fl=-CD*OZ’i t CKP+DX*$D 

Fl=Pl + (-OX->XP+r.Y*HP)*CO*CKP-DZ»XP*SD 
F1=F1+0X*CD* ( XN*CA-HN*SA ) +DZ*SD*( XN*CK-HN*$K> 
F2=-CDl>nZ4$KF-DY*SD 

F2=P2-*DY#CD*{XP*CKP+HP*SKP )-DZ*HP*SO 
F2=F2+OY*CD*{-XN*CA-*'HN*SA)+DZ*SP*(XN*SK+HN*CK) 
F3=0X* ( CP*CK*( -XN4CA+HN+SA )-SD*CK+XP*CD*CA) 

2 +DY4 ( -CP*'* K*(-XN>>C A+HN*9A ) +$D*SK— HP*CD*C A 1 

3 +DZ* l SD*( XP*CK+HP*SK-XN)+CC*CA ) 

F4=DX* tCD*SA* { XN*CK-HN*SK— XP) +SK*SD1 

2 +DY* (CD+SA* t-XN*SK-HN»CK+HP ) +CK*SD ) 

3 +PZ*( '0*(-XP*£K4HP*CK-HN)-CD#S A) 
FS=-{DX+SKP+DY*CKP)*CO 
F5=F5+PZ4C04 ( X P*SKP— HP*CKP ) 

F5=F5+DX*SD*(XN*SK +HN*CK l+DY^SD*! XN^CK— HN*SK ) 

F6 = TK*P5 
P7=CD*CKP 
F7=F7+XP*SP 

F7=F7+SP*(-XN*CK+HN*SK) 

F8=-CD*SKP 
Ffl=F8-HP*SD 
F6=F8+SD* ( XN*SK+HN*CK 1 
F9=SD 

F°=F,9+CD*(-XP*CKP-HP*SKP> 

F9=F9+CD*tXN>!‘C A-HN*$A 1 

F10=CD*(0X*SK.P+DY*CKP) 

F10=P10+DZ*CP*C-XP*SKP+HP*CRP) 

F10=F10+PZ*CDM-XN*SA-HN*CA) 

F11=(-DX*CKP+PY*SKP)*SD+PZ*CD 

F11=F 1 1 +C0* I DX*XP— DY4HP) ♦DZ*SD*(XP*CKP+HP*SKP) 

FI 1=F1 1 +DX4CD* (- XN#CK+FN*SK } +DY*CD* ( XN*SK+HN*CK ) 
2 -DZ*SD«MXN*CA-HN*XA) 

F 12=1 .DO 
F12=TK/24.D0 
C 

GlsDZ^CO^SKP 

G1 =G1 + DX*X p *CD*SKP-DY*H t, *CD*SKP 
G1=G1-DZ*SD*CXN*SK+HN*CK) 

G1 =OMG*C 1 
G2=-DZ*C0*CKP 

G7=G2+nY*Cn=*(-XP*SKP+HP*C t <P) 

C2=G2 +0 Z*SO* ( X N*CK— HN* SK ) 

G2 =n v iG*G2 

G3 = (UX*SK + nY-*CK)*SP 

G3=G3-{DX*SK+0Y>! ; CKl*Cn#(-XN*CA+HN*SA) 


21/55/22 
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0X60 

0161 

0162 

0163 

0166 

0165 

0166 
0167 
016 8 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 
017° 
01Q0 
0181 
0ie2 

' 0183 

0184 

0185 

0186 

0187 

0188 


0189 

0190 

0191 

0192 


G3=G3 + DZ«0*<-XP*SK+HP , 'I'CK> 

G3=0MG*G3 

G4=DX*CK*SD-DY*SK*SD 

G6=G4+DZ*SD*(-XP>!‘CK-HP*SK) 

G4=G4+DX*CD»-SA*(-XN*SK-HN*CK)+DY*CD*SA*{-XN*CK+HN*SK) 

G4=0MG*G4 

G5=-DX*C0*CKP4pY*Cn*SKP 

G5=G&40Z*C0*(XP*LKP4HP*SKP> 

G5=G5+OX*SD*tXN*CK-HN*SK>-DY*$D*(XN*SK4HN*CtO 

G5=!JMG*G5 

G6 =TK*G5 

G7=-CD*SKP 

G7=S7+S0*lXNvSK4HN*CK) 

G7=0MG*G? 

G8=-CD*CKP 

G8=GE* + SD*CXN*CK-HN*SK ) 

G8=DMG*G0 

G9=O.DO 

G9=G9+CD*(XP*SKP-HP*CKP) 

G°=0MG*G9 

G10=r>X*CD4CKP-DY*CD*SKP 
G10=G10~DZ*CP* (XP*CKP+HP*SKP> 

G10=UMG*G10 

G11=0X*SD*SKP+DY*SD*CKP 
G1 1=G1 1 +0Z*SD* t~XP*SKP+HP*CKP ) 

G1 1=01 1 4DX*CD * « XN* SK +HN*CK ) 4-D Y*C0* < XN*CK-HN*SK ) 

G1 1=DMG*G1 1 
G12=0.n0 
G13=l .00/24.00 
C 

C A. B ARE ROMS OF THE DESIGN MATRIX CORRESPONDING TO DELAY 

C AND DELAY DERIVATIVE OBSERVATIONS respectively 

C 

C INITIATE A AND B TO ZERO 

C 

DO 8 11=1. KK 

Min=o.oo 
Bfin=O.DO 
8 CONTINUE' 

C 

C SET PARTIAL S IN APPROPRIATE LOCATIONS OF A AND B ROWS 

C 


0193 


At 1)=F1 


019* 


A (2) = F2 


0195 


A {3) = F3 


0196 


. A ( 4) =F4 


0197 


At5)=F5 

o Q. 

0190 

C 

A < 6) = F6 

rd S 

019° 


&m=Gi 


0200 


B { 2 } =G2 

ip 

0201 


B (3t=G3 

0202 


9t4)=G4 


0203 


B ( 5 )=G5 

<3 to 

0 204 


B(6l=G6 

fe Q 


c 


cea 

0205 


11 = 3*144 


0 20 6 


12=1141 
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0207 


13=11+2 


0208 


A ( 11 ) =-F7 


o;o9 


A ( 12 )=— F6 


0210 

r 

A ( I3)=-F° 


0211 


B ( 11 ) =-G7 


0212 


B(I?)=-Ge 


0213 


B ( 13 1=-G9 


0214 


IC1=K+I 


0215 


IC2=1 Cl+IN 


0216 


A(IC1»=-F12 


0217 

r 

A( IC? )=-F13 


021? 

V 

B (IC1 ) =-Gl 2 


0219 

r 

8< IC2)=-G13 


C220 


J1=3*J+A 


0221 


J2=J 1+1 


0222 


J3=Jl+2 


0223 


A( J1 }= F7 


0224 


A(J?)= F8 


0225 


A ( J3 ) = F9 


0226 

c 

B ( J1 ) =G7 

o.g 

0227 


fi( 32 )=G8 


0228 

r 

B( J3)=G9 


0229 

V 

JCl=K+J 


0130 


JC2=JC1+IN 

v>* c> 

0231 


A ( JC1 )=F12 

& rd 

0232 

* 

A (JC2)=F13 


0233 

c 

B( JC1 )=G12 


0 234 

c 

8 ( JC2 )=G13 


0235 


IPI=3*IN+2*lP+5 


0236 


I P2=IP1+ 1 


0237 


A(IP1)=F10 


0238 

c 

A ( IP2 ) =F1 1 


0239 

V 

8 ( IP 1 ) =G10 


0240 

r 

B< I"2)=Gn 



v« 

c 

c 

FORMULATIOM OF 

L VECTOR ENTRY 

0241 

50 

0 SG=DX>) 1 F7+DY*F8+DZ*F9 

0242 


DSO=DSO +D0 t J ) -00 ( I) + ( DD t J ) -00 <!)>*( TK/24 .DO ) 

0243 

C 

L=CSO-OS 



C 

r 


C 

STORE PART 2ALS / NON-ZERO ELEMENTS OF H 8 ROWS ON DISK 8 FOR 


C 

LATER USE IN 

COMPUTATION OF RESIDUALS AND VTPV 


c 

c 

ADO CONTRIBUTIONS OF CURRENT OBSERVATIONS TO NORMAL EQUATIONS 


$**##*#*#**** + *********** 
C 
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0244 


WRITE (8, 351) F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11.F12,FI3,L 

0265 

351 

r 

FORMAT ( 10X*F40 .15 ) 

0246 

Lp 

DO 12 11=1, KK 

0 247 


U( 1I)=U(II )-A(II )*L 

0248 


DO 12 J J=1 ,KK 

0249 


NKEEP(II,JJ)=NKEEP(II,JJ )+ A( 1 1 ) *A { J J } 

0250 

12 

C 

CONTINUE 

0251 

c 

FRNG0=DX*G7+DY*G8+DZ*G9 

0252 


FRNGO=FRNGO+ (DO (J ) -DD ( I ) J/24.D0 

0253 

r 

L=FRNGO-FRNG 

0254 

L 

r 

WRITE (8,351 ) G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,L 

0255 

V 

DO 65 11=1, KK 

0256 


U(II)=UUI )-B(II)*l*PF 

0257 


DO 65 J J=1 ,KK 

0258 


N'KEEP ( 1 1 , J J )=NKEEP (II , JJ )+B( II)*B( JJ)*PF 

0259 

65 

f- 

CONTINUE 

0 260 

350 

f 

FORMAT (80X,2F19.9) 

0261 

W 

c 

GO TO 33 


C **** ** $*+ ****** **** **** X"***********************jJ:*#,)!**********::**$)*3} t :( 1 *#***:i:i)(***,)i* 

c 



c 

r 

SOLVE NORMAL EQUATIONS, COMPUTE PARAMETER CORRECTIONS 


V 

C**************************************************** ************* He************* 
C 

0262 

c 

66 

DO 1 1=1,8 

0263 


DO 1 J= 1 ,KK 

0264 


IK=I+KK 

026 5 


NKEEP ( I K, J )=CC (I , J ) 

0266 


NKEFP (J,IK)=NKEFP(IK,J) 

0267 

1 

r 

CONTINUE 

0268 

L 

r 

DO 18? IW=1,I0PT 

0269 

L 

DO 183 1=1, K8 

0270 


DO 1P3 J=! ,K8 

0271 


N( I, J)=NKEEPCI ,J) 

0272 

183 

r 

• CONTINUE 

0273 

V 

IN1=7+IN*3 

0274 


00 2 I=IN1,K 

0275 


N(I,I )=N( I, I )+PRATIO(IW) 

0276 

2 

r 

CONTINUE 

0277 

V 

r 

CALL DMINV(N,Kfi,DET,Ll,L2) 

0278 

V 

00 187 1=1, KK 

0279 


XX ( I ) =0 . DO 

0280 

18? 

C 

CONTINUE 
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0261 
0 282 
0283 
028 A 


0285 

0286 

0287 

0288 

0289 

0290 

0291 


02 ° 2 
0293 

02« A 
0295 

0 296 

0297 

0298 

0299 
0 300 

0201 

0302 

0303 

0304 
0 305 
0206 
0307 
0208 
0-09 
0210 


0311 

0312 

0313 

C31 A 
0*15 
0316 


PO 14 1=1, KK 
PO 14 J=1,KK 
XX <1 )=XX(I )+N( I,J)*U< J) 
1A CONTINUE 

C 

C WRITE XX 

C 

00 9 1=1,5 
CALL REOPI(XXd)) 

9 CONTINUE 

C 

ISTART=7+3*IN 

C 

00 26 I=ISTART,K 
CALL REOPKXXd)) 

26 CONTINUE 

C 

C 


C 

C READ A AND B { ROWS OF DESIGN MATRIX) FROM DISK 8, 

C FORMULATE RESIDUALS AND COMPUTE VARIANCE COVARIANCE MATRIX 

C 

c 

85 REWIND 8 

VTPV=O.ODO 

C 

VABS1 =0 .DO 
VA6S2=0 .DO 


C 

C 


C 


23 

C 


DO 22 II=l,ICOUNT 
READ ( 8, 100 ) 1,J,IP 

©>© 



DO 27 IDUM=1,2 

READt 8,351) ( F < J J ) , JJ =1 , 13 ) ,L 
V=0.D0 

Htf® 

§1 

DO 23 J J=1 , 6 
V=V+F (JJ)*XX(JJ) 

CONTINUE 

feS 



11=2*1+4 
12=11 + 1 
13=11+2 

^55 


Jl=3*J+4 


J 2=J1 +1 


J3=J1 +? 

V = V+F(7)*{XXCJ1>-XXCU)) +F(B)*(XX<J2)-XX<I2) > 
2 +F(9)<=(XXU3)-XX(I3)1 


IP l=3*IN+2*IP+5 


IP?=IP1+1 

V=V+F { 10 )*XX {IP1)+F(11)*XX(IP2) 

IC1=K+I 
IC 2=IC 1+TN 
JC1 =K+ J 
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0317 


JC2=JC1+IN 

0318 

r 

V=V+F(12)*(XX( JCl)-XX(ICl) )+FU3)*<XX(JC2)-XX(IC2)l 

0319 

c 

V=V+L 

0320 


I F t I0UM.EG.2 > GO TO 28 

0221 


VABS1=VARS1+DABS(V) 

0222 


VTPV=VTPV+V**2 

0323 


GO TO 27 

0324 

u 

28 

VABS2=VABS2+DABS<V) 

0325 


VTPV=VTPV+PF4V**2 

0326 

27 

f* 

CONTINUE 

0327 

27 

r 

CONTINUE 

0328 


VABS1 =VABS1/IC0UNT 

0329 


VAtsS2=VABS2/IC0UNT 

0330 


ITN=3*IN+6 

0331 

r 

IIN1=IIN+1 


v» 

c 

ALTERNATIVE FOR SO COMPUTATION 


c 

SO=VTPV+DFLOAT ( IM*2)*( SDBD/R100 )*92 

■ 

c 

f* 

•SO=SO/DFLOATCICOUNT*2-KK+2*IM) 

0332 


DO 172 I=ITM1,K 

0333 


VTPV=VTPV+PRATI0(IW)*XX{I)**2 

0334 

•172 

CONTINUE 

0 335 

r 

S0=VTPV/DFL0AT(2*IC0UNT-KK+2*IH) 

0336 

V 

DO 24 1=1, KK 

0337 


DO 24 J=1,<<K 

0338 


N( I»J)=NtIfJ)*SO 

0339 

24 

C 

CONTINUE ' 

0340 

C 

Kl=K+l 

0341 


KIN=K+IN 

0342 


KIN1=KIN+1 

0343 


DO 41 J=1,KK 

0244 


N(1,J)=R100*N(1,J> 

0345 


N(2» J )=R100+N( 2»J ) 

0346 


N{ 3, J)=RI00*N{ 3, J) 

0347 


N(4,J}=RIOO*N(4,J) 

0248 


N«5, J)=RASM3*N<5,J> 

0349 

r 

N(6, J)=RHR0M9*N(6,J) 

0350 


DO 42 1=7, IIN 

0351 


N < I , J J =N{ T , J ) *R1 00 

0352 

42 

r 

CONTINUE 

0353 

U 

DO 43 1= I INI »K 

0354 


N(I»J)=NtI»J)*RA$M3 

0355 

43 

CONTINUE 

0356 

' 

DO 47 1=K1,KIN 

0357 


N(I,J)=(N(I,J) /CLI GHT 1 *1 000.00 

0358 

47 

CONTINUE 
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0359 

0560 

0361 

0362 

0363 

0264 

0265 

0366 

0367 
0566 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 
0281 

0382 

0383 

0384 
0285 

0386 

0387 

0388 

0389 

0390 
0291 
0392 


0393 

0394 


0295 

0296 


C 


48 

C 

41 

C 


C 


45 

C 


46 

C 


49 

c 


52 

C 

44 

C 

227 

229 

222 

C 


86 

C 


696 

C 

8 50 


C 

356 


00 48 I=KIN1,KK 

/CL1GHT J *1000.00 

CONTINUE 

CONTINUE 


CO 44 1=1, KK 

N(I,1)=R100*NCI,1) 

NtI,?)=R100=*NII,2) 

NII,3)=R100*N(I,3> 

NCI,4)=RIOO*N(I,4) 

N<Z,5)=RASM3*N(I,5) 

N( I,6)=RHRDM9*N(1,6) 

DO 45 J=7, IIN 
N{ I,J)=N(I,J)*R100 
CONTINUE 

DO 46 U=IIN1,K 

N( I,J)=N(I,J)*RASM3 

CONTINUE 


«>© 

q 

l! 

3sS 


00 49 J=K1 »KIN 

MI,J) = IN(I,J) /CLIGHT) *1000 .DO 
CONTINUE 

DO 52 J=KIN1,KK 

N{ T,J)=(N( I, J ) /C LIGHT ) *1000. DO 

CONTINUE 


CONTINUE 


FORMAT (10X»3F20.4) 
FORMAT ( 10X»2F22.6) 
F0RMAT(5X,F40.9) 

on 86 1=1, KK 

N( I , 1 )*D5QRT (N (I , I ) ) 

CONTINUE 


ICNT2=ICPUNT*2 
I0F=ICNT2-KK+2*IM 
WRITE (6,696) I CNT2 ,IDF 

FORMAT (/////5X,»N0 OF OBSERVATIONS =*,I5//5X 
2, ‘DEGREES CF FRF EDOM =‘,I5/) 

WRITE (6, 850) I EXP , SDBD,SOBF , SDQUA ( I W ) 

FORMAT { 1 l‘///15X, ' EXPERIMENT NO =',I5///10X 

1, 'A PRIORI STANDARD DEVIATIONS * ,//5X 

2, ‘DELAY DISTANCES (CM ) * , 10X, F 15 .5/5X 

3, ‘DISTANCE DERIVATIVES ( M/HOUR )‘, IX, F15.5/5X 

4, ‘QUASAR COORDINATES ( CM ) » ,7X ,F15 .5/ ) 

WRITE ( 6 ,356) { N( I , I) , 1=1 ,6 ) 

FORMAT {///20X, ‘STANDARD DEVI ATI ONS */20X 
1, • A******************* //5X 

?r • POLA R MOTION = • ,?X, F14 . 5 , 3X t * CM • ,/26X ,FI4.5r 3X, ‘CM V5X 
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2, ‘PRECESSION-NUTATION s • ,F 14. 5.3X , • CM • ,/26X .F14.5.3X. * CM* , 

4, * TH F T A ZERO = • , 9X , F14. 5 ,3X, • ARCSEC/1000 »/5X 

5, »0MG = , ,16X»Fl4.5,3X f ‘REV PER DA Y/10**9 >//10X 

6, 'STATION COORDINATES (X f Y,Z) IN CM‘/> 

0397 


WRITE 16 »357 ) (MI.I >,I-7,IIN> 

0396 

357 

FORMATf 10X,3F?0.5) 

0399 


WRITE (6 »35P 1 

0400 

358 

FORMAT { /10X. "QUASAR COORDINATES (RA,DECI IN ARCSEC/IOOO*/) 

0401 


WRITF( 6.359) ( N( I , I ) , I =1 INI . K ) 

0402 

3 59 
c 

FORMAT ( 10X.2F20.5 > 

0403 

W 

WRITE ( 6.361 > 

0404 

361 

FORMAT (//10X t * STATION CLOCK OFFSETS IN MSEC*/) 

0405 


WRITE (6 .362 ) (NfI,I),I=KI»KIN) 

0406 

362 

F0RMAT(?X,F?0.9) 

0407 


WRITE ( 6 >363) 

0408 

363 

F0KMAT(//10X, • STATION CLOCK DRIFTS IN MSEC PER DAY' /) 

0409 

C 

WRITE (6.362 i t Nf I, I )» I=KIN1»KK) 

0410 

C 

182 

C 

86 ’ 

CONTINUE 

C411 

STOP 

0412 


END 



3. Name of Program: PROG02 
Function: 

Adjustment of observations with a time span equal to the station step. 

Station coordinates are considered constant over the whole set of observa- 
tions, while earth rotation parameters are treated as constant but 
different over a prearranged set of subintervals. 

Input: 

Standard deviations of observations and radio source coordinates, approximate 
values of parameters, station and radio source identification numbers, 
distance and distance rate observations. 

Output: 

Standard deviations of adjusted parameters. 

Subroutines required: 

MADD, MATPRO, TRAPRO, PRO TR A for matrix operations and RED PI . 

(Following is a listing of PROG02 and supporting subroutines, except for 
REDPI already included in VLBI SIMULATOR.) 
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FORTRAN 


0001 

0002 


0003 


XV G1 RELEASE 2.C MAIN DATE = 76300 21/56/32 

I,******************************************************************************* 




********** 

c********* 

PR0G02 

********** 

c********* 


********** 

c********* 

VLB I ADJUSTMENT PROGRAM 

********** 

c* *v****** 


********** 

C********* 

THIS VERSION PERFORMS THE ADJUSTMENT OF 

********** 

c********* 

A SET OF OBSERVATIONS TREATING STATION 

********** 

C#***v**jI <* 

COORDINATES AS CONSTANTS. EARTH ROTATION 

********** 

c********* 

PARAMETERS ARE CONSTANT BUT DIFFERENT 

********** 

c********* 

OVER SUBINTERVALS OF THE ORIGINAL TOTAL 

********** 

c********* 

TIME INTERVAL OF THE AVAILABLE OBSERVATIONS. 

********** 

C********* 


********** 

c********* 

OBSERVATIONS: DELAY DISTANCES AND 

********** 

C 

THEIR DERIVATIVES 

********** 

C**#* ** v** 


********** 

t ********* 

UNKNOWNS: EARTH MOTION PARAMETERS 

********** 

c ********* 

STATION C QUASAR COORDINATES 

********** 

C********* 

CLOCK OFFSETS AND DRIFTS 

********** 

c ********* 


********** 

(,********* 

SYSTEMS DEFINITION: 1 

********** 

c********* 

OUASAR FIXED: WEIGHTED QUASAR COORDINATES 

********** 


EARTH FIXED: BY MEANS OF INNER CONSTRAINTS 

********** 

c < ******** 


********** 

t********* 


********** 

c ******* ** 

SUBROUTINES REQUIRED : 

********** 

c ********* 

REDPI 

********** 

c******* ** 

MADD 

********** 

C******** * 

MATPRO 

********** 

C ********* 

TRAPRO 

********** 

C*** |<***** 

PRQTRA 

********** 

£********* 

. 

********** 


C ****** ********************* *********** ***************************************** 
O******** 1 ********** 

IMPLICIT REAL*? (A-H»L-Z) , 

DIMENSION DO (3 ),DD(3) , X t 3 ) ,Y ( 3 ) ,Z <3 ) » RA ( 9) , 0 (9 ) ,ED( 33,8) ,U (39) 

2, X'KEEPl 'i°,ZQ ) , A(2«) ,3(3« ) , NB {6,33), UDI C 33 ) , HI ( 6, 33 ) ,RI (33,33) 

3,0H«1,M) ,KLI <33),H{6,33) ,RL<33) ,UD < 33 > ,HTEDO ( 33 , 8 ) , UK (41 ) 

4»L?K ( 41 ) ,0,1' <33 ,32 ) , FB ( 33 , 8 ) , GBL <33 ) ,FET<33) ,GU< 33 ) , FBE0T( 33,6) 

5, C-6HT (33,6) ,PPTN< 33,6) ,XB I < 6 , 20 ) , STD < 153 ) ,XB (33 ) 

6, TO ( 6 , 8 ) , UDDI to) , NOD < 6, 6 ) , LI 6( 6) ,L26(6),TI (6) ,RK<6,6) ,RT<6) 

7»RKE00(6*8) ,ETKE(F,8) ,FF{8,8) »ET0T(8),FSDT(8),FBTL(’fi)»F8TU(R) 

f » SDL 1 6 ) ,FDF(6,f) ,EFET<6,6) ,GDO< 6, 6 ) ,HGHT( 6, 6) ,HFEN (6,6 ) , EFEN < 6,6 ) 
9,NFFFN(6,6),HIX8(6 ) ,NSOL (6),F (13) [ 

C DIMENSION DO(IN),DD(IN),X(IN),Y(IN) , Z ( IN ) ,R A( IM,D ( IM) , ED (K ,8 ) 

C 2,U(KK> ,NKECP < K K, KK ) , A ( KK ) , 8 ( K.K ) »NB( 6* K)UDI(K) , HI ( 6, K ) , RI (K ,K ) 

C 2,GE(K8,KP) ,RLI (K),H(6,K) ,RL < K > , UD (K ) ,HTEOO ( K,8 ) ,L IK (K 1 ,L2kfK> 

C 4,GP<K,K),FCMK,6),G3L(K),FET(K),GU(K),F8E0 < K ,6 > ,GBHT (K ,6 ) 

C £,FFTN(K,6) ,XS I (6 , 1 DCNT ) • STD (6*IDCNT +K ) , X8 <K ) » 

C 6,rD0( 6,fi),l!DDI (6), NOD (6, 6) ,L16(6) ,L26(6) ,TI (6),RK(6,6) ,RT(6) 

C 7,RKFOO(6»8),ETKE(fl»8} , FF ( 8 ,8 ) ,ET0T { 3 ) , FEDT( 8 ) , FBTL ( 8 ) , FBTU (8 ) 

C 8, SOL 16) ,r.DF(6,€),EFET(6,6) ,GDD(6,6) ,HGHT ( 6, 6 ) , HFEN ( 6 ,6 ) , EJf EN ( 6,6 ) 

C 9 ,N£Fr.N (6,6)HIX8(6) ,NSfiL(6) ,F( 13) 

C WHERE: K=5*IN+2*IM,KK=K+6,KB=K+8,IN=NO OF STATIONS, IM=NO OF QUASARS 

EQUIVALENCE (HT FDO, F3 ) , ( UDI , LIK, GBL > , [ RLI ,L2K, FET ) , ( RI ,G8 ) 

2, ( RKEDO , EDF,GBHT) , ( EDO ,FETN) , ( RT,HIXB) , (L16,TI ) , ( L26,NSQL> 

3, ( RK,NEFSN>, (HGHT, HFEN, EFEN), <ETKE,FF), (ETOT, FBTL, FBTU) 



FORTRAN 


0004 

0005 

0006 


0007 

0008 

OC 09 
0 010 


0011 

0012 


0013 

0014 

0015 

0016 
0017 
0 0 IF 
GO 19 
0020 
0021 
0022 

0 02 3 
0024 


0025 

0026 
0027 
002 ? 

0029 

0030 
0C31 


0032 
0 033 
0024 
0035 


IV G 1 RELEASE 2.0 MAIN DATE = 76300 21 / 56/32 

C REAP EXPtKTMENT IDENTIFICATION NO 

RFAD { 5 , 100 ) IEXP 
WRIT ? ( 6 * 477 ) I EXP 

977 FORMAT ( • 1 1 //// 10 X , ‘EXPERIMENT NO =*,I 5 / 10 X 

?, • ********<*********** t , ///// ) 

C 

c 

C RFAD A PRIORI STANDARD DEVI ATIONS t 

C SOSO : OF DELAY DISTANCE IN CM 

C SDB F : OF DELAY DISTANCE DERIVATIVE IN M/HOUR 

C SOQUA: OF QUASAR COORDINATES IN CM 

C 

RE AD { 5 r 101 ) SDBDtSDBF 
101 FORMAT ( 5F1 6. 5 ) 

C 

WMTE< 6 , 240 } SDbD, SDBF 

2-0 FORMAT! 10 X, 'OBSERVATIONAL NO I SE « / 15 X, * D I STANCES 
2 , C 10 . 5 / 16 X, 'DISTANCE RATES ',F 10 . 5 //) 

C 

RF AD ( 5 » 1 01 ) SOQUA 
PRATI 0 =SPeP** 2 /SDQLA **2 
C 

C READ EARTH STEP INTERVAL IN HOURS 

C 

RF AD ( 5 * 101 ) FRSTEP 
/RITE ( 6 , 243 ) ERSTEP 

243 FORMAT (//IPX, 'EARTH STEP INTERVAL (HRS) =*,F 10 . 2 /> 

C 

R = 637 1000 . DO 
C L IGHT = 3 .DO* ( 1 0 . 00 ** 8 ) 

ClIGHT=CLICMT/R 

PI= 4 .D 0 *L'ATAN( 1 .D 0 ) 

R 10 G=R* 100 .DO 

KHKPM 9 = ( 12 .D 0 /PI )* ( 10 . 007 * 9 ) 

RASM 3 =( 180 . DC»/PI >* 3600000 . DO 
C 

PF = ( S 0 ! 5 D/RlOO)/(SD 8 F/R ) 
p F=PF **2 

C PF IS NOW IN (EARTH RADII/HGUR)**- 2 * PRATIO IN (EARTH RADII)**-? 

C 

C RF AT IN= NO OF STATIONS C IM= NO OF QUASARS 

£**** *.!$!***$***#>**#******:************************* ******************************* 

RFAO ( 5 . 1 O 0 ) IN, IM 
ICO F ORMAT ( 1 C 15 ) 

WRITE ( 6 , 920 ) I N, IM 

920 FORMAT ( 5 X, • NO OF STATIONS =*,I 5 / 5 X,»N 0 OF QUASARS =*,I 5 ) 

K= 5 *IN+ 2 *IM 

Y K=K + 6 
KF=K +8 

c 

C 

DO 31 1 = 1 , IN 
00 ( 1 ) = 0.00 
00 ( 11 * 0.00 
51 CONTINUE 
C 

•C READ APPROXIMATE VALUES OF PARAMETFRS THETA ZERO, OMG, 

(;^<**ii< >!i*r,*X<*** ************** ***************************************************** 
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0036 

0037 

0038 
0C39 
00 AO 

0041 

0042 


CC43 
0044 
004 5 
0046 
0 047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


0056 

C057 

0058 

0059 

0060 
0061 

0 062 

0063 

0064 


0065 

0066 

0067 

0068 
006 9 
0070 
C071 

0072 
0 073 

0074 

0075 

0076 


C xm, Y(I), Z(I), (1 = 1, IN), RA ( J ) , D(J), <J=1,IM) 

C ANGLES ARi IN DFGREES, DISTANCES IN METERS AND OMG IN RADIANS PER HOUR 

Rf AD (5,101) TH 
WRIT E ( 6 ,F60 ) TH 

660 FORMAT ( //5X* 'THETA ZERO = *,F20.9/) 

TH£T=TH*PT/180.DO 

0MG=PI/12.00 

C 

WR ITF ( 6 » 450) 

450 FORMAT <////, 10X,* STATION C RADIO SOURCES APPROXIMATE COORDINATES* 

?,//3X, • STATIONS', 10X, 'LONG* 

2, 12X, *LAT* ,11X, 'R* ,/) 

C 

DO 5 1=1, IN 

RF AD (5, 701) I DM, I LI, I L2, ALONG, APHI,AR 
WRITfc( 6,701) I, XL 1, I L?, ALONG t A PHI ,AR 
701 FORMAT! IX, T2, IX, A4,A2,3F15.I> 

APHI=APHI*PI/1PC.D0 
AL0NG=AL0NG*PI/1B0.D0 
C CONVERT TO CARTESIAN 

Xd )=AR*DCOS(APHI)*DCOS( ALONG) 

Yd)=AR*DCOS(APHI)*DSIN( ALONG) 

Zd ) = AR*DSIN(APHI) 

C CONVERT TO EARTH RADII UNITS 

X(I)=X(I)/R 
Y(I)=Y(I)/R 
Z( I)=Z ( I )/R 

5 CONTINUE 
C 

C 

WRITF (6,451) 

451 F0RMAT(//11X, 'RADIO SOURCE 10X ,» RA *, 17X , 'DEC » ,//) 

C 

DO 6 J=1,IM 

READ! 5,321) 1DM , 1 L X , I L 2, IL3 , RA{ J ) ,D ( J ) 

WRITE (6, 321) J,IL1,IL2,IL3,RA(J),D(J) 

321 FORMAT (1X,I3,7X,3A4,F17.7,E20.9) 

C CONVERT TO RADIAN UNITS 

RA (J)=RA (J )*PI/160.D0 
DU) = D(J)*PI/180.D0 

6 CONTINUE 
C 

C 

C FORMULATE CONSTRAINT MATRICES ED, EDO 

DO 11 J = 1,S 
DO 9 1=1,6 
EDOC I , J 1 =0 .DO 
9 CONTINUE 

DO 11 1=1, K 
ED (I , J ) =0.00 
11 CONTINUE 

C 

F DO (1,51=1. DO 
£00(2 ,4 ) =1 .DO 
EDC (5,6 )=— 1 .DO 
C 

DO 3 1 = 1, IN 
I e 1=3* 1-2 
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0 C 77 

ocas 

0079 

ooeo 

ooei 

0052 

0053 
OOB A 
OOP 5 
0 086 
008 7 

or.p*' 

0069 

0090 

0091 

0092 
0 C93 
0094 


0095 
OC 96 
0097 
Ov.98 
Ot)99 
0100 

0101 

0102 

0103 

010s 

C 2 05 
0106 

0107 

0108 
0109 


0110 

0111 

0112 

0113 

0114 

0115 


I £2=1 21 + 1 
1S3=I62+1 
ED(I°1,1 > =1 .do 
£D(IE-?, 21 = 1.00 
£01153,31=1.00 
ED 1152 ,4 >=— Z 11) 

£0(163,41=7(1) 

EDUS1 ,5) = zm 
80(193 , 5 )=-X (T 1 
E0(:31,6)=-Y(I) 

E0(IB?,6)=xm 

ICl=3*IU+2*Tf' + I 
1C2=IC 1+IN 
ED ( IC1 , 7 1=1 .DO 
EDI IC2 , 81=1.00 

3 CO., T INU£ 

ICCNT=0 

Icnu.MT=o 

C 

C INITIATE H, RK, RL, RT , UD, OB TO ZERO 

C +** + ** :p:t* + ****V**W* + * + *+************+******>i<***********:)::S‘*’4‘’S‘***l)‘***>)‘*l)‘ #*#****#* 

C 



DO 71 1=1, K 




KL (I 1=0 .00 




UD ( 1 1 =0 .DO 




CO 71 J=l,6 

O 

g 


Hi J, I 1=0. DO 

hrj 

3 

71 

CONTINUE 

►d 

'O 

C 

DO 72 1=1,6 

81 


RT (I }=0 -DO 

P3 

H 


DO 72 J=1 ,c 

fO 

k 

T) 


RK ( I , J 1 =0 .(iC 

Cj 

7? 

CONTINUE 

> 

P> 


00 72 1=1, KB 
DO 73 J=1,K6 

1 

s 

ra 


OH I, J) =0.00 

kJ 

}— i 

CO 

73 

CONTINUE 




c 

c 

C INITIATE NKEFP.U TO ZERO 

C+*** *+***+>)** ***+****+***)<.)‘***>!‘>S<****>)‘**if<**+*<****’«‘****i)'*’)'***>)f***+'M‘*'M'-****>)t+’)' 

55 DO 16 1=1, KK 

11(11=0.00 
DO 16 J=1,KK 
NK E EP ( I » J 1 =0 .00 
If CONTINUE 

C 

TH=TH£T+ (EkST£P/l2.D0)*PI*0FL0AT( IDCNT ) 

C 

C + + + + + ++ + + *+++*+:), ++*!+«+ *+++ + *++*++*3)$ 

S* 

C READ OBSERVATIONS AND IDENTIFIERS: 

L 0S= DISTANCE IN DETERS OF IJ BASELINE OBSERVATION TO QUASAR IP 

C AT TINE TK (IN HOURS] AFTER SOME INITIAL EPOCH TO 

C++++* * +*+ + ,<*++* + !»**#+* + *+W + ++-)r:t::'ii+>(i+i!t:i< + +++++ + **+ + :(<+i):i)::{c* + :(c:Mii)ti(t*=jt+i(i+!*c:)t>((+>)ti)i!4ci|iijci(ti)t*>!<i+ 

L 

f 
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0116 

0117 


0118 

0119 

0120 
0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 
0123 

0134 

0135 

0136 


0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0184 

0155 

0156 

0157 

0158 

0159 

0160 
0161 
0162 

0163 

0164 

0165 


C 

33 

104 

C 


O O 

. 

•n.8 

gg 

l§~ 

3s 


READ{5,104> It J,IP,DS»FRNG r TK,ICHEK 
FORMAT < 3 15, 3F1 5. 5, 215) 

ICHEK IS A CODE INDICATING END OF DATA FOR ICHEK=1 
IF (ICHEK .6Q.1) GO TO 66 
IF (ICHEK .E0.2 ) GO TO 66 
I COUNT = I COUNT ♦ 1 
TK=TK-FRSTEP*DPLOAT (IDCNT ) 

OS=t)S/R 
FRNG=FRNG/R 
DX=X( J)-X(I) 

PY=Y( J)-Y(I) 

DZ = Z(J)-2(n 

CA=DCOS (RA (IP ) ) 

SA=DSIN(RA(IP) ) 

C0=nC0S (DIIPJ) 

SD=0$IN(P< IP) ) 

Y1=0M6*TK+TH 

Y2=Y1-RA(IP) 

CK=PCOS ( Yl ) 

SK=[)$ IN ( Y1 ) 

CKP=PC0S(Y2) 

SKP=DAIN(Y2) 

FORMULATE PARTIALS 

Fla-CD^PZ^CKP+DX+SD 

F 2 ==“C D Z =f SK P-GY^ S 0 

F3=-DX*SD*CK+DY*SD‘»SK+DZ*CD*CA 

F4=DX*SK+SD+DY*CK*SD-0Z*CD*SA 

F5=-(PX*SKP+DY*CKP )$CD 

F6= TK*F5 

F7*CP+CKP 

FC=-CO*SKP 

C 9=SD 

F10=CD4(DX*SKP+DY*CKP) 

F11={-0X*CKP+DY*SKP)*SD+DZ*CD 
F 1 2- 1 . DO 
F 1 3=TK/24 .00 

G1=DZ*CP*SKP 
Gl=nr,G«Gl 
G2=-DZ*CC’*CKP 
G2 apMG*R2 

G3=(DX*£K+DY*CK)*SD 

G3=OMG*G3 

G4 =DX * C K * S D-0 Y*S K * S 0 
G4a0MG*G4 

G5 a-OX^CD^CKP+DY^CP^SKP 

G5=0f'G*G5 

G6=.TK*G5 

G7=-CD*SKP 

G7=0MG*G7 

GP=— CD’l'CKP 

GP =OMG*GP 

G9=O.DO 
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0166 


G10=DX*CC*CK P— OY*CD*SKP 

0167 


Gl 0=0MG*G1 0 

0 1 66 


Gll=OX*Sn*SKP+rjY*SD*CKP 

0169 


GU-OMG^Gll 

0170 


G12=0.D0 

0171 

z' 

G13=l .D0/24. DO 


c 

r 

INITIATE A AND B TO ZERO 


C ********** S*******#***#****#****##***********#*****************^************** 

r 

0172 


DO 8 II=1»KK 

0173’ 


At II)=O.DO 

0174 


ptin=o."DO 

0175 , 

8 

CONTINUE 


C 

c 

SET PARTI AL S IN APPROPRIATE LOCATIONS Or A ROW 


C******* ****** ********** tf**************************************************** 1 *** 
n 

0176 

V 

A ( 1 )= FI 

0177 


A (? ) =F2 

0178 


A ( 3) -F3 

0179 


A(4)=94 

oieo 


At 6) = F5 

0161 

r 

A 16) =F6 

0182 

U 

B 111 = 01 

0183 


B(2)“G? 

0134 


8 ( 3)=G3 

0185 


e (4}=64 

oieo 


8 (5}=G5 

0187 

r 

R ( 6 ) =G6 

0188 

V 

11=3*1+4 

0189 


12=11+1 

0190 


13=11+2 

0191 


A ( 11 ) =-F7 

0192 


A ( I? ) =-FP 

0193 

c 

A t 13 )=-F9 

0194 


B( Il) = -G7 

0195 


B ( 12 ) =-G8 

0196 


B ( 13 ) =-G9 

0197 

L 

IC1=6+3*IN+?*IM+I 

01 Q 9 


1C?=IC1+1N 

0199 


At IC1 >=-F12 

0 2C0 

r 

A(IC2)=-F13 

0201 

L> 

8 ( IC1 ) =-Gl 2 

0202 

r 

Et(IC2)=-G13 

0203 


J 1 =3* J+4 

0204 


J2 =J 1+1 

0205 


J3=J1 +2 

0206 


A ( Jl ) = F7 

0207 


A l J2 )= F 8 

0208 


A { J3 ) = F9 


C 
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0209 

0210 
0211 

0212 

0213 

0214 
02 IS 

0216 

0217 

0218 

0219 

0220 
0221 

0222 

0223 


0224 

0225 

0226 


0227 

0228 

0229 

0230 

0231 

0232 

0233 
0 234 


0235 

0236 

0237 


023e 

0239 

0240 

0241 

0242 

0243 


0244 


B( J1)=G7 
S t J2 )=G0 
B( J3)=G9 
C 

JC1=6+3*IN+2*IM+J 
JC2=JCI+IN 
A ( JC 1 ) =F 32 
At JC2)=FI3 
C 

6 ( JC 1 ) = G12 
6 { JC2 ) =G 13 
C 

IP 1=3*1 N+2*I P+5 
I P2=IP 1+ 1 
AtlPl >=F10 
A ( IP2)=F1I 
C 

B ( I P 1 ) =% 1 0 
B t IP2 ) =G1 1 
C 

C PORMULATIOM OF t VECTOR ENTRY 

C 

50 DS0=DX*F7+DY*F8+DZ*F9 

DSO=OSO+DO(J)-DOtI J+{DD( J J-DD(I ))*lTK/24 .DO) 

L=DSO-DS 

C LW 1=L*R100 

C 

c 

C £$)!■ *$ + « * 4 # ** **>*>>: ***** * + ifit **$ + ****>> A* !>******$* + ******+*<< ******* 

c STORF A, B > L AND OBSERVATION IDENTIFIERS ON DISK 8 FOR LATER 

C USE IN COMPUTATION OF RESIDUALS AN COVARIANCES 

C ****** ********&##>»**$*****##* **>S'**#*##*************##***+#****:4‘ !&*************♦* 

c 

WRirE( 8 , 10 A) ICHEK,I,J,IP 

WRITE ( 8,351 ) F1,F2,F3,F4,F5,F6 t F 7,F8,F9,F10,F11,F12,F13,L 
351 F0RMAT(10XtF4O .15) 

C 

rn 12 IT=l,hK 
u(n)=u(in+Atii)*L 
DO 12 JJ=1,KK 

NKEEPt II,JJ)=NKEEP(II,JJ ) + At 1 1 )*A ( J J ) 

12 CONTINUE 

C 

C 

FRNG0=DX*G7+DY*G£+DZ*G9 

F R NGO= F RNGO * ( D D < J ) ~D D ( I ) ) / 24 . DO ??????? 

L=FRNGO-FRNG 
C LW 2=L*R 

C WRITE ( 6 * 350) LW1.LW2 

C 

WRITE (8,351) GltG2tG3»G4»G5»G6*G7»G8TG9TG10»Gll»Gl?.G13»L’ 

DO 65 11=1, KK 
U{II)=UtII)+P(II)*L*PF 
DO 64 JJ=1,KK 

NKEEPtll, JJ)=NKEEP til , J J ) + B ( 1 1) *B fJJ)*PF 
65 CONTINUE 

C 

350 FORMAT (60X»?F19.9) 
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0245 


0246 

0247 

0248 

0249 
025 0 

0251 

0252 

0253 

0254 
02 5 5 
02*6 


0 2 57 

0258 

0259 
0 260 


0261 
0262 
026 3 

0264 

0265 
0 266 
0267 
0269 
0269 


027 C 

0271 

0272 
0272 
C274 

0275 

0276 


0277 

0278 

0279 

0280 

02ei 
0282 
0 283 
0284 


C 

GO TO 33 
C 

c 

C 

C 

66 WRIT E ( $ » 100 ) ICHEK 
C 

DO 20 1=1,6 
U0DI{1)=U( I) 

DO 20 J = 1 » 6 
NDPU ,J)=NKESP(I,J) 

20 CONTINUE 

DO 21 1=1,6 
DO 21 J=1,K 
j6=J+6 

NBU, J->-=NKfcEPt I, J6 > 

21 CONTINUE 
C 


DO 2? 1=1, K 
16=1+6 

UOKI )=U(I6J 

22 CONTINUE 
C 

CALL DMINV tNC'D ,6,DFT, L16 ,L26 > 

CALL MATPRO C NDD»NB ,HI,6,6,K) 

CALL TRAPR0(NB,HI,RI,K,6,K1 
DO 23 1=1, K, 

DO 23 J=1,k' 

16=1+6 
J6=J+ 6 

QB ( I , J ) =QP- 1 1 , J ) +NK EEP ( 16 , J61-R1 ( I , J ) 

23 CONTINUE 
C 

CALL MATPRC(NDO,UODI,TI t 6,6,1) 

CALL TRAPR0(N5,TI,RLI,K,6,1> 

CALL MA D0( F, HI ,6 ,K ) 

CALL MADDtRK , NOD, 6, 6) 

CAlL MAOD t RL.RLI *K ,1) 

CALL MADI t RT ,T I, 6 , 1 ) 

CALL MAC'D ( UD, UDI, K, 1 > 

C 



C STORE UDDI ,NDD ,HI ,TI ON DISK 9 

C** *<* *.*>!■ *#*>>* *#****##=X*#4**!»*#:fc***#4**>»*****#***l*iM!##*!jl*#!«' #****# ********«:******* 

C 

DO 25 1=1,6 
DO 2 C J=l,6 
WRITE (9,351) NDC ( I , J ) 

25 CONTINUE 
C 

DO 26 1=1,6 
00 26 J=1,K 
Wk 3TE ( ° ,35 1 ) HI (I , J) 

26 CONTINUE 
L 


0285 


DO 27 1=1,6 
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0186 


WRITF (9*351) Till) 

0287 

27 

CONTINUE 

0288 


IDCNT=IDCNT+1 

0239 


IFnCHEK.EQ.2) GO TO 67 

0290 

C 

GO TO 55 


C 

c 

0291 

c 

67 

r 

CALL TRAPROtH, EDO, HTEOO t K ,6 t 8 ) 

0292 

%# 

DO 28 1=1, K 

0293 


00 28 J= 1 , 8 

02° A 


JK=J+i< 

0295 


QS(I,JK)=ED(I» J l-HTEDO ( I , J ) 

0296 


0*(JK,I)=QP<I,JK) 

0297 

28 

c 

CONTINUE 

0296 

V 

CALL MATPROCRK, EDO »RK EDO, 6,6,8) 

0299 

r 

LALL TRAPRO(EDO,RKEDO,ETKE,8,6,8) 

0300 

V 

00 29 1=1,8 

0301 


1K=T*K 

0302 


00 29 J= 1, B 

0303 


JK=U+K 

03C A 


0E>tIK,JK)=-E7KE{I,J) 

0305 

29 

r 

CONTINUE 

0 306 

V 

1 PF=3# IN+1 

C307 


IPL=3*IN'+2*rM- 

0308 


DO 3A I=IPF,IPL 

0309 


QB(I,I)=CB(I,1 J+PRATIO 

0210 

34 

r 

CONTINUE 

0311 

V 

r 

CALL DMI\VtOB,K8,DET,LlK,LZK) 

0312 

u 

00 30 1=1, K 

0313 


00 30 J=1,K 

031A 


GP ( X » J )=CB 1 1 , J ) 

0315 

30 

r 

CONTINUE 

0216 

U 

DO 31 1=1, K 

0217 


PO 31 J=l,8 

0318 


JK=J+K 

0319 


FB (I » J }=Q3 1 1 , JK ) 

0320 

31 

c 

CONTINUE 

0221 

L 

00 32 1=1,8 

0322 


IK=I+K 

0323 


CO 32 J=1 , 8 

0224 


JK=J+K 

0J25 


FFII , J 1=08 { IK , JK ) 

0326 

32 

r 

CONTINUE 

0327 

u 

CAlL MATPRO{GB,RL,GBL,K,K,l) 

0328 


CALL TRAPRO(EOO,RT,ETOT,0,6,1) 

0329 


CALL MATPR0(FS,ET0T,FET,K,8,1) 
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0330 

r 

CALL MATPR0(CB,UD,GU,K,K,1 ) 

0331 


DO *1 1=1, K 

0332 


XB (I )=G£L( 1)+FET(I)-GU(I) 

0333 

*i 

r 

CONTINUE 

033* 

v» 

CALL M*TPR01PF,ET0T,FEDT»8,8,1) 

0335 


CALL TRAPROf FB,RL,F6TL»8,K,11 

0336 


CALL NAPP(F60T,FBTL,8,1) 

0 337 

r 

CALL TRAPR0(FB,UD,FBTU,8,K,1) 

0336. 

i* 

no ** 1 = 1 , e 

0339 


FEPT ( I)=PEPTII ) -FBTU( I ) 

03 AO 

44 

r 

CONTINUE 

03*1 

L 

CALL MATPRO ( EPO , FEDT ,SOL ,6,8,1) 

03*2 


CALL PRQTR A( FB ,ED0 , FBEOT ,K ,8 , 6) 

03* 3 


CAuL I / ATPRn(cD0,FF,£DF,6,8,6) 

03** 

r 

CALL PROTRA(EOF, E0O,EFET ,6,8,6) 

0 3*5 

c 

REWIND 9 

03*6 

V 

c 

DO 60 IPAY=1,IPCNT 

0347 

c 

DO 35 1=1,6 

03*8 


CO 35 J =1 , 6 

03*9 


RE An t ° , 351 ) NDD(I , J) 

0350 

35 

r 

CONTINUE 

0351 

L 

DO 36 1=3,6 

0352 


DC 36 J= 1 * K 

0253 


RF AO < 9 , 3 51 ) HICI.JJ 

035* 

36 

r 

CONTINUE 

0355 


DO 37 1=1,6 

0356 


* EAD ( 9,351) TI(I) 

0357 

37 

r 

CONTINUF 

0356 

U 

DO 38 1=1,6 

0359 


00 38 J=l,6 

0360 


GDPt I , J )=NPD (I ,J) 

0361 

38 

r 

CONTINUE 

0362 

u 

CALL PROTR A t GB »HI , GBHT ,K, K ,6 ) 

0 363 


CALL NATPRO(HI ,GBHT, HGHT , 6 ,K , 6 ) 

036* 


CALL MAPP ( GO!r, HGHT ,6,6) 

036 5 


CALL MAT C R0(FBE0T,NDD,FETN,K,6,6) 

0266 

r 

CALL MATPkO(HI ,FETN,HFEN,6,K,6) 

0 367 


DO 39 1 = 1, '6 

0368 


DC1 39 J=l,6 

0 369 


GDD( I , J )=GDD ( I , J (+HFEN ( I, J )+HFEN( 

0370 

3* 

CONTINUT 

0371 


CALL MATPR0(6FET»NDD»EFEN»6,6»6) 

0372 


CALL MATPR0(NDD,EF£N,N£FEN»6,6*6J 

0373 

. 

CALL MAPD{GDD,NEFEN,6,6) 


C 
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0374 

0375 

0376 

0377 


0370 
0 3 79 


03P0 

0301 

0382 

0303 


0384 

0365 

0386 

0367 


0383 

03S9 

0^90 


0391 

0?92 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 

0403 

0404 

0405 

0406 

0407 

0408 

0409 

0410 

0411 

0412 

0413 

0414 

0415 


C 

c 

c 

Of 48 1=1,6 

IPF=6*IDAY-6+I 

STO{TPF)=GDD(I»I) 

48 CONTINUE 
C 

C 

CALL MAT PRO { NDP, SOL , NSOL ,6,6, 1 ) 

CALL MATPR0(HI,XB,HIXB,6,K,1) 

C 

00 40 1=1,6 

XSI(I,IDAY )=**T 1(1) ~NSQL( I ) -HI XB ( I ) 

40 CONTINUE 

60 CONTINUE 
C 
C 

C s^*******##* *$*#**.* j********** #****:jc* #* *$$*** *************** ********** ********** 

C 

C 

00 49 1=1, K I 

IPL=6*IPCNT+I 

STD(IPL)=GB(I,I) 

49 CONTINUE 
C 

REWIND 8 
ID AY=1 
VTPV=0.D0 

c 

77 READ (8,100) ICHEK,I,J,IP 
IF { ICHCK .F Q.l ) GO TO 88 
I F ( TCH £K . SO. 2 ) GO TO 89 
I S 1=3*1—? 

IS2=IS1+1 

I S3=I S2+1 

JSl=3*J-2 

JS2=JS1+1 

JS3=JS2+1 

IP1=3*IN+?*IP-1 

IP?=IP1+1 

ID 1 = 341 N+24IM + I 

J01=3*IN+2*IM+J 

ID2=I01+IN 

JD2=JD1+IN 

DO 45 JJ=1,2 

RE AD ( 8 , 35 1 ) ( F ( II ) ,1 1 =1 , 13 ) , L 
V = L 

DO 42 11=1,6 

V=V+F (II J’XXO I t II ,IDAY) 

42 CONTIMUF 

C 

V=V+F (7>«(XB(JS1>-XB(IS1 ) ) +F ( 8 )* ( XB { J S2 )-XB { IS2 ) ) 

2 + F ( 9 ) *< XB ( J S3 )-XB ( I S3 ) ) 

V=V+F(10)*XB(IP1)+F(11 )*XB { I P2) 

V=v+F(12 )4(X8( JDl)-XB(IDin+F(13>*( XB ( JD2 ) — XB ( ID2 ) ) 

IF (JJ.E0.2 ) GO TO 43 
C VW=V*R100 
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C 

WRITE (6,453) VW 

0416 


VTPV =VT PV+V**2 

C417 


GO TO 45 

0418 

43 

VTFV=VTPV+PF*V**2 


C 

VW=V*R _ 

WR ITE { 6 »454) VW Q. S 

FORMAT (50X »F30 .15) ^ p. 


C 


C454 


C453 

FORMAT 12 OX r F 30. 15) *1> S 

0419 

45 

CONTINUE Q2, - 

0420 


GO TO 77 S ^ 

0421 

86 

IDAY=I0AY+1 ^ CT 

0422 / 


00 TO 77 

0423 

89 

IPF=2*TN*1 C3 

C424 


IPL=3*IN*2*IM Q 

0425 


DO 46 I=IPF,IPL C© 

0426 


VTPV=VTPV+PRATI0*X6(I )**2 *3 t* 
CONTINUE 4 s * u 

0427 

46 

r 

0428 


IS?=2*IC0UNT 

0420 


JS2=5*IN+6*I0CNT 

0430 


I S3=IS2-JS2 

0431 


1DC=(?4.D0/ERST£P+0.1D0) 

0432 


IOC =1 DC NT/ IOC 

0433 


WRITE ( 6 1 242 ) I PC ,1 52 , JS2 , IS3 

0434 

24? 

FORMAT (///10X, 'DAYS OF OBSERVATION^ ,I5/10X, 'NO OF OBSERVATIONS = » 
2,I5/10X,'N0 OF UNKNOWNS = », 15/1 OX, • DEGREES OF FREEDOM **♦15//) 

0435 

r 

S0=VTPV/DFL0AT(IS3) 

0436 

u 

IPL=K+6*I0CNT 

0437 


DO 47 1=1, IPU 

0438 


IF(STD( D.GT.O.DO) GO TO 999 

0439 


WRITE ( 6,996 ) I.STD(I) ' 

0 44 0 

998 

FORMAT ( /////•***** ERROR **44* I =* , I 5/10X • STD= * ,F90 .30/////) 

0441 

099 

STOCI >=DS<3RT(STD(I )*S0) 

0442 

h7 

r 

CONTINUE 


V 

‘ C********* WRITE STANDARD DEVIATIOS 

0443 


WR ITE ( 6,888 ) 

0444 

888 

r 

FORMAT ( *) *//10X» 'STANDARD DEVIATIONS */10X, » ****************♦**»// 
2,3X, 'DAY'tIOX, 'POLAR MOTION' , 13X, 'PRECESSION - NUTATION', 7X 

3 , 'THE T A ZERD * , 7X , ' OMEGA' /12X, *KSI { CM ) • ,7X, *ETA (CM)',7X 

4, *KSI ( Cr 1 } » , 7X , • ETA ( CM )', 9X ,•( CM )* , 6 X, • (REV/DA Y) *10**9 »//) 

0445 


DO 85 IDAY=1»IDCNT 

0446 


I PL= 6 *I PA 7 

0447 


IPF=IPL-5 

0448 


DO 86 11 = 1,5 

0^49 


JJ=IPF+II-1 

0450 


CALL REDPI (STD (JJ) ) 

0451 


STDUJ )=STD( JJ)*R100 

04 52 

86 

CONTINUE 

0453 


STO(IPL)=STD(I PL) *RHRDM9 

0454 


WRITE (6,889) IDAY, (STD ( J ) , J=IPF, I PL ) 

0 455 

889 

FORMAT (/1X,I4,6F15 .5) 

0456 

85 

CONTINUE 

04 47 


WR ITF ( 6 ,221 ) 

0458 

221 

FORMAT ( ' l'///10X, 'STATION COORDINATES • »//,17X, , X',9X,*Y' ,9X 


2, *Z‘/15X, 1 (CM) *,6X,' (CM) *,6X, '(CM)'//) 
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0459 


IPF=6*IDCNT+1 

OAc'C 


IPL=6*I0CNT+3*IN+2#IM 

0461 


II=IPL-2*IM+1 

0462 


DO 91 I=II,IPL 

0 463 


CALL RFDPI (STD ( I ) ) 

0464 

91 

C 

CONTINUE 

04 t 5 

DO 75 I=IPF,IPL 

0466 


STD!I)=S"ID{I)*R10Q 

0467 

75 

r 

CONTINUE 

046$ 

L 

IPL=IPL-2*1M 

0469 


WRITE (6,222) ( STD (I) , 1=1 PP , IPL) 

0 470 

222 

c 

FORMAT { 10X,3F10.2) 

047 1 


WRITE ( 6 t 273 ) 

0472 

223 

FORMAT ( ///// 10X , *OUAS AR COORD INATES * ,//l 6X, *RA « ,8X , »DEC • ,/15X 
2t • (CM)*,6X,*(CM) •//) 

0473 


IPf=IPL+l 

0474 


IPL=IPL+2*IM 

0473 


WRITE (6, 224 > < STD(I) » 1=1 PF , IPL) 

0476 

224 

r 

FORMAT ( I0X»2Flo.2 ) 

G*t77 


WRITE" ( 6 » 225 ) 

0478 

?25 

c 

FORMAT (/////10X, ‘STATION CLOCK OFFSET? AND DRIFTS V13X 
2, • (MSEC) »,6X, • (MSEC/DAY) •//) 

0479 

L 

DO 76 1=1, IN 

O-^RO 


I C 1= I PL + 1 

0461 


IC?=IC1+IN 

rp4£2 


STDdCl ) = ( STD ( IC1 )/CL I6HT )*1000»D0 

0463 


STDEIC2)=(STr { IC2)/CLIGHTm0CO.OO 

04 64 


WR ITC (6,226) STD ( IC1 ) , STD (IC2 ) 

0*»65 

276 

FORMAT (4X,2F16. 12) 

0486 

76 

c 

CONTINUE 

0467 

L 

STOP 

0466 


DEBUG SUBCHK 

04 6 9 


END 
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0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 
0009 


SUBROUTINE MADDCA,B,L,M) 

#*■ 

Kj 

C 

C A = A + B 

C L#M L*M L*M 

C 
C 

IMPLICIT REAl*8(A-M,Q-Z) 
DIMENSION A(L,M) ,B (L»MJ 
DO 5 1=1, L 
DO 5 J=1,M 
A(I,J)=A(I,JJ+B(I,JJ 
5 CONTINUE 
RETURN 
END 
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or.oi 


0002 

0003 

0004 

0005 

0006 

0007 

0008 
000 ° 
0010 
0011 


SUBROUTINE XATPROI A, B »C t L ,M, N ) 
C 
C 

C C = A # B 

C L*N l*M M*N 

C 

c 

IMPLICIT RF4L*8(A-H,0-Z) 
DIMENSION A ( L»M) »B(M,N)»C(L,N) 
DO 5 1=1 , L 
DO 5 J=I,N 

cti) j)=o.no 

DO 5 K=1,M 

C(I,J)=C(I,J)+A(I,K)*B(K,J) 

5 CONTINUE 
RETURN 

Fkc 
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0001 


0002 
0003 
0 004 
0005 
OCOto 
0007 
OOOP 

0009 

0010 
OOU 


SUBROUTINE TRAPROtA»B t C*L*M»N) 
C 
C 

C T 

C C = A * B 

C L*N L*M M*N 

C 
C 

IMPLICIT REAL*e(A-H, 0 -Z) 
DIMENSION A(M,U,B(M,N>,C{L»N) 
DO 5 1 = 1 , L 
DO 5 J= 1 ,N 
C(I, J)=O.DO 
DO 5 K= 1 ,M 

C{I,J)=CtI ,J)+A(K,I)*B(KtJ) 

5 CONTINUE 
RETURN 
END 
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i 


C001 


0002 

0003 

0004 

0005 
OOOA 
COO? 
0008 
000 ° 
OCIO 
0011 


SUEROUTI N£ PROTRAI A,B ,C,L,M T N) 
C 
C 

C T 

C C = A # B 

C L*N L*M M*N 

C 
C 

IMPLICIT RFAL*8(A-H,0~Z) 

DIMFMSION A(L,M)*6«N,M),C(LtN) 

00 5 Iwl,L 

00 5 J = 1,N 

CtI,J)=0.D0 

00 5 K«=1,M 

C< I,J)=C(I» J)+A<I,K)*B(JtK) 

5 CONTINUE 
RETURN 
SNC 



